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CONVERSION  TABLE 


Conversion  factors  for  U.S.  Customary  to  metric  (SI)  units  of  measurement 

MULTIPLY  BY  TO  GET 

TO  GET  BY  DIVIDE 


angstrom 

1.000  000  X  E  -10 

meters  (m) 

atmosphere  (normal) 

1.013  25  X  E  +2 

kilo  pascal  (kPa) 

bar 

1.000  000  X  E  +  2 

kilo  pascal  (kPa) 

barn 

1.000  000  X  E  -28 

meter2  (m2) 

British  thermal  unit  (thermochemical) 

1.054  350  X  E  +  3 

joule  (J) 

calorie  (thermochemical | 

4.184  000 

joule  (J) 

cal  (thermochemical )/cm5 

4.184  000  X  E  -2 

mega  joule/m  (MJ/m  ) 

curie 

3.700  000  X  E  +1 

*giga  becquerel  (GBq) 

degree  (angle) 

1.745  329  X  E  -2 

radian  (rad) 

degree  Fahrenheit 

t,  =  (t*F  +  459 . 67 ) / 1 . 8 

degree  kelvin  (K) 

electron  volt 

1.602  19  X  E  -19 

joule  (J) 

erg 

1.000  000  X  E  -7 

joule  (J) 

erg/second 

1.000  000  X  E  -7 

watt  (W) 
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3.048  000  X  E  -1 
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foot-pound-force 
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gallon  (U.S.  liquid) 

3.785  412  X  E  -3 

meter  (m  ) 

inch 

2.540  000  X  E  -2 

meter  (m) 

jerk 

1.000  000  X  E  +9 

joule  (J) 

joule/kilogram  (J/kg)  (radiation 

dose  absorbed) 

1.000  000 

Gray  (Gy) 

ki lotons 

4.183 

terajoules 

kip  (1000  lbf) 

4.448  222  X  E  +3 

newton  (N) 

kip/inch2  (ksi) 

6.894  757  X  E  +3 

kilo  pascal  (kPa) 

ktap 

newton-second/m2 

1.000  000  X  E  +2 

(N-s/m2) 

micron 

1.000  000  X  E  -6 

meter  (m) 

mi  1 

2.540  000  X  E  -5 

meter  (m) 

mile  (international) 

1.609  344  X  E  +3 

meter  (m) 

ounce 

2.834  952  X  E  -2 

kilogram  (kg) 

pound-force  (lbs  avoirdupois) 

4.448  222 

newton  (N) 

pound-force  inch 

1.129  848  X  E  -1 

newton-meter  (N*m) 

pound-force/inch 

1.751  268  X  E  +2 

newton/meter  (N/m) 

pound-force/foot2 

4.788  026  X  E  -2 

kilo  pascal  (kPa) 

pound-force/ inch2  (psi) 

6.894  757 

kilo  pascal  (kPa) 

pound-mass  (lbm  avoirdupois) 

4.535  924  X  E  -1 

kilogram  (kg) 

pound-mass-foot2  (moment  of  inertia) 

ki logram-meter2 

4.214  011  X  E  -2 

(kg*m2) 

pound-mass/foot* 

kilogram/meter 

1.601  846  X  E  +1 

(kg/m1 ) 

rad  (radiation  dose  absorbed) 

1.000  000  X  E  -2 

“Gray  (Gy) 

roentgen 

coulomb/ki logram 

2.579  760  X  E  -4 

(C/kg) 

shake 

1.000  000  X  E  -8 

second  (s) 

slug 

1.459  390  X  E  +1 

kilogram  (kg) 

torr  (mm  Hg,  0“  C) 

1.333  220  X  E  -1 

kilo  pascal  (kPa) 

*The  becquerel  (Bq)  is  the  SI  unit  of  radioactivity;  1  Bq  =  1  event/s. 
“The  Gray  (Gy)  is  the  SI  unit  of  absorbed  radiation. 
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SECTION  1 
INTRODUCTION 


Trans  ionospher ic  communications  systems  which  rely  on  the  propa¬ 
gation  of  radio  frequency  (RF)  signals  may  suffer  severe  performance 
deqradation  when  the  ionosphere  is  disturbed  by  high  altitude  nuclear 
explosions.  PRPSIM  (Properties  of  Radio  Wave  Propagation  in  a  Structured 
kmized  Medium)  is  a  FORTRAN  77  computer  program  which  has  been  designed 
to  compute  a  set  of  parameters  to  enable  the  user  to  assess  the  severity 
of  propagation  effects  for  arbitrary  link  conf igurat ions  in  simulated 
nuclear  environments.  PRPSIM  has  been  developed  under  the  sponsorship  of 
the  Defense  Nuclear  Agency  (DNA),  and  it  uses  the  official  DNA  propagation 
algorithms  to  compute  the  signal  parameters.  These  algorithms  have  been 

derived  in  a  number  of  references.  The  purpose  of  this  volume  (Volume  II) 
is  to  describe  in  detail,  and  in  one  source,  all  the  computations  per¬ 
formed  by  PRPSIM.  For  instruction  in  running  the  PRPSIM  code,  the  reader 
should  refer  to  the  user's  manual  (Volume  I). 

The  severity  of  the  disruption  of  communications  systems  in 
nuclear  environments  depends  on  many  factors,  but  these  factors  may  be 
grouped  into  two  general  categories:  (1)  phenomenology  parameters  such  as 
the  number  of  bursts,  weapon  yields,  altitudes,  and  magnetic  locations  of 
the  bursts,  and  (2)  link  parameters  such  as  signal  carrier  frequency, 
antenna  character i st ics ,  and  link  geometry  and  velocity  relative  to  the 
nuclear  environment.  PRPSIM  is  a  system  propagation  effects  code  which 
has  been  designed  to  provide  the  user  with  maximum  flexibility  in  the 
specification  of  link  parameters.  It  has  no  capability  for  providing 

phenomeno loqy  parameters,  however.  Instead,  PRPSIM  has  been  designed  to 

operate  in  simulated  nuclear  environments  generated  by  the  MRC 
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phenomenology  code  SCENARIO.  SCENARIO  is  a  physics  based  code  whose  simu¬ 
lations  of  disturbed  environments  are  obtained  as  the  solutions  of 
coupled,  time-dependent  equations  which  represent  the  dominant  physical 
and  chemical  processes  operating  in  the  disturbed  ionosphere.  SCENARIO 
was  developed  under  DNA  sponsorship  and  was  desiqned  specifically  to  model 
nuclear  environments  generated  by  multiple,  high  altitude  explosions  above 
100  km.  The  operation  of  the  SCENARIO  code  is  beyond  the  scope  of  this 
documentation,  but  detailed  descriptions  are  available  in  several  of  the 
references  listed  at  the  end  of  this  volume.  Phenomenology  data  files  are 
written  by  SCENARIO  at  discrete  times  after  the  initial  explosion.  PRPSIM 
may  then  interface  with  the  SCENARIO  environment  at  these  specific  times, 
or  at  any  intermediate  time  by  interpolating  the  phenomenology  data 
between  bracketing  times. 

The  radiant  energy  generated  by  nuclear  explosions  produces  very 
highly  ionized  regions  in  the  atmosphere.  The  magnitude,  location, 
spatial  extent,  and  duration  of  the  ionization  are  all  sensitive  functions 
of  the  weapon  yield,  altitude,  and  other  phenomenology  parameters.  At 
relatively  low  altitudes  (below  about  50  km),  most  of  the  weapon  energy  is 
initially  deposited  in  a  relatively  small  reqion  around  the  explosion, 
producing  a  compact  but  intensely  ionized  "fireball".  At  higher  altitudes 
(above  100  km),  the  weapon  enerqy  is  dispersed  over  a  larger  volume,  pro¬ 
ducing  a  loss  intense  but  much  larger  and  lonqer- 1  ast  inq  region  of  ioniza¬ 
tion.  Another  phenomenon  which  occurs  for  high  altitude  explosions  and 
which  has  an  adverse  effect  on  RF  signal  propagation  is  the  formation  of 
structured  ionization  in  the  plasma.  Typically,  the  plasma  structure 
consists  of  lonq  filaments  of  high  electron  density  aligned  along  the 
geomagnetic  field  lines.  These  filaments  (or  "  str i at  ions" )  may  extend  to 
v°r y  high  altitudes  and  produce  a  scintillation  of  RF  signals  analogous  to 
the  scintillation  of  light  waves  propagating  through  a  thermally  turbulent 
atmosphere.  Atmospheric  disturbances  typical  of  both  low-altitude  and 
h  iqh-a  1 1  i  tude  explosions  are  illustrated  in  Figure  1.  The  user  should  be 
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Figure  1.  Propagation  through  a  striated  region. 


aware  that  PRPSIM  does  not  account  for  many  of  the  effects  typical  of 
low-altitude  and  surface  explosions. 

The  resulting  signal  degradation  may  be  grouped  into  two  classes 
of  effects:  (1)  mean  effects  such  as  absorption,  phase  shift,  time  delay, 
etc.  that  are  caused  by  the  enhanced  mean  ionization,  and  (2)  stochastic 
effects  such  as  amplitude  and  phase  scintillation,  arrival  angle  scatter, 
and  time  delay  jitter  which  are  caused  by  the  plasma  structure  and  its 
motion  relative  to  the  propagation  path.  The  various  types  of  propagation 
disturbances  and  their  overall  effects  are  summarized  in  Table  1.  The 
significance  of  each  of  these  propagation  effects  on  the  performance  of  a 
particular  system  depends  strongly  on  the  link  design  (carrier  frequency, 
data  rate,  type  of  modulation  and  coding,  tracking  loop  configuration, 
etc.)  and  on  the  power  margin  and  system  requirements.  Some  degradation 
of  the  system  performance,  relative  to  that  in  an  undisturbed  environment, 
is  inevitable  in  a  nuclear  disturbed  environment.  Careful  attention  to 
the  nuclear  propagation  threat  and  application  of  practical  engineering 
techniques  can  provide  hardened  link  designs  with  acceptable  levels  of 
survivability  in  most  nuclear  environments.  The  purpose  of  PRPSIM  is  to 
provide  the  system  engineer  reliable  data  with  which  to  assess  system 
survivability  in  specific  nuclear  environments. 
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TABLE  1.  Summary  of  various  types  of  propagation  disturbances  and 
their  overall  effects. 


PROPAGATION; 

DISTURBANCE 

EFFECTS  ON  PROPAGATED  SIGNAL 

EFFECTS  ON  SYSTEM  OPERATION 

Mean  Ionization  Propagation  Disturbances  calculated  by  PRPS1M. 

Abscrpt  ion 

Loss  of  RF  energy  due  to  electron- 
ion  and  electron-neutral  colli¬ 
sions  in  ionized  regions,  causing 
signal  attenuat ion . 

Reduction  of  signal -to-noise 
ratio,  causing  possible  loss  of 
data  and  disruptions  of  signal 
acquisition  and  tracking. 

No  i  se 

Random  RF  energy  emission  from 
hot,  ionized  fireball  regions. 

Increased  background  noise, 
causing  reduction  of  siqnal-to- 
noise  ratio  and  possible  loss 
of  da. a. 

Phase  and 

Dodo  1 er 

Shift 

Signal  phase  variation  due  to 
changes  in  TEC  alona  propagation 
oath  . 

Possible  loss  of  phase/frequency 
lock  and  difficulty  in  signal 
acquisition  and  tracking. 

Ti^e  De'ay 

Signal  time-of -arrival  variation 
due  to  changes  in  TEC  along  pro¬ 
pagation  path. 

Possible  loss  of  timing  lock  and 
difficulty  in  bit  or  PN  code 
acauisition  and  tracking. 

Faraday 

Rotat  ion 

Sianal  polarization  rotation  due 
to  maqnetoionic  anisotropy. 

Severe  ar n p  1  itude  fading  with 
linear  polarization;  no  effect 
with  circular  po 1 ar izat ion . 

Striated  Ionization  Propagation  Disturbances  Calculated  by  PRPSIM. 

Arp!  itude 
Scintillation 

Random  fluctuations  of  signal 
amplitude  due  to  scattering  of 
signal  energy  by  ionization 
irregularity  structure. 

Severe  amplitude  fading,  causing 
increased  data  error  rates  and 
possible  disruption  of  signal 
acquisition  and  tracking. 

Phase 

Scintillation 

Random  fluctuations  of  signal 
phase  due  to  TEC  variation  and 
scattering  of  signal  energy  by 
ionization  irregularity  structure. 

Loss  of  signal  coherence,  caus- 
possible  loss  of  phase  lock  and 
catastrophic  failure  of  demodu¬ 
lation  . 

Anqij  1  ar 
Scatter  i no 

Random  fluctuations  of  signal 
ang  1  e-of -am i va  1  due  to  scatter¬ 
ing  of  signal  enerqy  by  ioniza¬ 
tion  irregularity  structure. 

Reduction  of  signal -to-noise 
ratio  due  to  scattering  loss 
in  large  antennas,  causing 
possible  loss  of  data. 

Frequency 

Se  1  ec  t  i  ve 

Sc int i 1 1  at  ion 

Random  distortion  of  signal  wave¬ 
form  and  signal  time  delay  due  to 
frequency  dependent  angular 
scatter  ing . 

Degradation  of  signal  detection, 
tracking,  and  demodulation; 
possible  loss  data  due  to  inter¬ 
symbol  interface. 
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SECTION  2 

RF  SIGNAL  PROPAGATION 


2.1  INTRODUCTION  TO  THE  MODELS  OF  PROPAGATION  DISTURBANCES. 

The  primary  task  in  any  analysis  of  the  propagation  of  an  RF 
siqnal  through  a  turbulent,  ionized  medium  is  a  quantitative  description 
of  the  variation  of  the  ionization  (or  equivalently,  of  the  index  of 
refraction)  as  a  function  of  time  and  location  within  the  medium.  All 
PRPSIM  calculations  are  based  on  the  assumption  that  the  fluctuations  at 
any  particular  time  may  be  described  by  a  spectral  density  function  in 
three  dimensions  which  decreases  with  increasing  wave  number  according  to 
an  inverse  power  law.  The  power  spectral  density  is  parameterized  by  data 
generated  by  the  SCENARIO  environment  simulation  code  at  each  point  in  the 
SCENARIO  grid.  The  relevant  parameters  are  transformed  to  a  coordinate 
system  which  is  invariant  with  respect  to  location  within  the  grid,  and 
the  set  of  quantities  which  describe  the  scintillation  of  the  RF  signal 
are  then  calculated  usinq  the  equations  derived  by  Wittwer  in  the  reports 
DNA-5304D  1  and  DNA-IR-82-02 2. 

PRPSIM  is  designed  to  operate  in  an  anisotropic  propagation 
environment.  The  scale  lengths  of  the  fluctuations  in  three  directions 
are  required:  parallel  to  the  magnetic  field,  and  in  two  directions 
perpendicul ar  to  the  field.  The  scale  lengths  are  presumed  to  vary  with 
location  in  the  grid.  The  spectral  index  of  the  power  law  can  have  any 
value  between  1.0  and  6.0.  Currently,  the  value  of  the  index  is  assumed 
to  be  invariant  with  respect  to  location,  but  in  the  future  the  code  will 
be  generalized  to  accommodate  a  variable  index. 
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2.2 


SCINTILLATION,  ABSORPTION,  AND  MEAN  ELECTRON  DENSITY  EFFECTS. 


The  set  of  parameters  required  to  describe  the  propagation  of  a 
scintillated  RF  signal  through  a  turbulent  ionized  medium  is  derived  and 
discussed  by  Wittwer  in  DNA-5304D 1  and  DNA-IR-82-02 2.  For  completeness, 
we  summarize  here  the  equations  used  by  PRPSIM  to  calculate  these 
parameters.  The  reference  numbers  of  the  equations  in  Wittwer' s  documents 
are  provided  also,  should  the  user  have  any  questions  regarding  the 
derivations.  All  parameters  are  calculated  by  the  module  RFPROP.  All  the 
data  required  to  define  the  propagation  path  and  its  end  members  (the 
transmitter  and  receiver)  are  passed  directly  to  RFPROP  by  PRPSIM  as 
arguments.  These  data  consist  of  the  transmitter  and  receiver  location 
and  velocity  vectors,  the  number  of  carrier  frequencies  at  which  simula¬ 
tions  are  desired,  the  array  containing  those  frequencies,  and  an  array 
containing  the  effective  antenna  beamwidtbs  at  those  frequencies. 


The  derivation  of  the  parameters  is  based  on  the  assumption  that 
the  fluctuation  of  the  electron  density  and  the  index  of  refraction  along 
the  path  can  be  represented  by  a  power  spectral  density  (PSD)  function 
whose  form  is  a  simple  power  law: 


PSD(Kr,Ks,Kt)  = 


8  .i  V 2  ^ 2  L rL  SL t  r(n)/r(n-3/2) 

(1  +  l2K2  +  l2*2  +  L,2-K2)n 
r  r  s  s  t  t 


if  (K  K2+  K2)  i2  <  1 


(2-1) 


otherwise 


where  K r  and  are  the  spatial  wave  numbers  of  the  fluctuations  in  two 
orthogonal  directions  perpendicul ar  to  the  geomagnetic  field  vector,  and 
is  the  wave  number  in  the  direction  parallel  to  the  field;  An?  is  the 
variance  of  the  index  of  refraction  about  its  mean  value;  r(n)  is  the 
gamma  function  of  argument  n;  and  n  is  proportional  to  the  spectral  index 
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( 2-2) 


n  =  —  spectra]  index  +  1 
2 

(See  Equation  1  of  DNA-IR-82-022.)  This  form  of  the  PSD  is  consistent 
with  the  observation  of  very  sharply  defined  large  scale  structure  in  the 
ionosphere.  The  spectral  index  (or  n)  can  be  interpreted  as  a  measure  of 
the  sharpness  of  the  structure.  We  have  assumed  that  the  spectral  index 
of  the  turbulence  does  not  vary  along  the  path.  We  have  fixed  its  value 
at  2  because  in  situ  measurements  and  theoretical  analyses  imply  that  this 
value  best  character izes  the  local  structure  of  the  ionsophere.  (See,  for 
example,  References  1  through  10  of  DNA-IR-82-022.)  The  spectral  index 
may  be  assigned  any  value  within  the  range  1.0  to  6.0  (1.5  <  n  _<  4.0), 
however  .* 


Lr,  l_s,  and  Lt  are  the  striation  outer  scale  sizes  as  measured 
perpendicular  and  parallel  to  the  field  vector,  and  l  is  the  inner  stri¬ 
ation  scale  size.  (These  parameters  are  assigned  by  the  data  retrieval 
submodule  DATAPT.)  In  the  context  of  the  PSD  (Equation  2-1)  used  through¬ 
out  this  documentation,  the  scale  size  L  is  the  distance  (measured  in  a 
particular  direction)  over  which  the  argument  of  the  PSD  (KL)  changes  by  1 
radian  in  that  direction.  In  general,  the  outer  scale  size  L  is  not  iso¬ 
tropic  in  direction  because  the  large  scale  turbulence  is  strongly  depen¬ 
dent  on  the  geomagnetic  field.  The  inner  scale  size  i  is  the  smallest 
distance  over  which  fluctuations  can  be  detected.  Its  value  is  usually 
considered  to  be  isotropic  in  direction  because  the  small  scale  turbulence 
fluctuates  rapidly  and  is  not  strongly  influenced  by  the  geomagnetic 
field.  The  structure  of  a  plasma  in  a  magnetic  field  is  characterized  by 
long  striations  aligned  along  the  field  lines.  Because  the  striations  are 


Wittwer's  derivations  are  valid  for  the  range  1.2  <  n  <  4.0.  The 
differential  phase  variance  (equations  2-37  and  2-38)  is  undefined  for 
n  £  1.5,  however.  Until  a  valid  expression  for  the  phase  variance  is 
derived  for  n  £  1.5,  we  shall  restrict  the  value  of  n  to  be  greater 
than  1.5. 


long  compared  to  their  transverse  dimensions,  OATAPT  assigns  the  longi¬ 
tudinal  scale  size  l_t  to  a  constant  value  of  150  km  for  all  space 
throughout  the  seguence  of  the  simulation.  The  character ist ic  transverse 
scale  sizes  are  generated  by  the  SCENARIO  code.  Typically,  these  are  of 
the  order  of  10  km.  OATAPT  assigns  the  (isotropic)  inner  scale  size  a 
constant  value  of  10  meters.* 

All  the  parameters  required  to  parameterize  the  scintillated  RF 
siqnal  are  functions  of  integrals  which  are  integrated  over  the  length  of 
the  path  from  the  transmitter  to  the  receiver.  RFPROP  utilizes  a  simple 
integration  technique  to  evaluate  each  of  these  integrals  --  with  one 
exception  which  will  be  discussed  in  detail  later.  First,  a  set  of  path 
segments  (a z -j )  is  defined;  second,  a  set  of  data  describing  the  environ¬ 
ment  at  the  midpoint  of  the  seqment  is  retrieved  for  each  of  the  segments; 
finally,  the  integrals  are  summed  with  the  assumption  that  the  value  of 
the  integrand  at  the  center  of  the  segment  is  a  good  estimate  of  the  mean 
value  of  the  integrand  over  the  length  of  the  segment;  i.e.,  we  compute 
the  integral  of  each  integrand  f  as 

zend  N-l  z-,,+z-  N 

r  f(z)dz  -  Y  f(_l±li)(z.+1-zi)  *  Y  f(<z.>)Az.  (2-3) 
o  i=0  2  1  i=l  1 

where  N  is  the  number  of  segments  along  the  path.  When  we  represent  the 
integrands  by  functions  which  can  be  integrated  analytically,  it  can  be 
shown  that  the  accuracy  of  this  "mean-rectangular"  technique  is  superior 
to  the  trapezoidal  technique,  and  that  it  is  comparable  in  accuracy  to 
Simpson's  technique. 


The  values  used  for  the  longitudinal  outer  and  the  isotropic  inner 
scale  sizes  were  suggested  in  a  private  communication  with  L.  A. 
Wittwer  of  PNA. 
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The  set  of  path  segments  and  the  points  at  which  the  integrands 
are  calculated  are  generated  by  the  submodule  PATHPT  (described  in  Section 
0 . .  The  set  of  points  and  line  segments  is  such  that  each  point  corre¬ 
sponds  to  t he  midpoint  of  each  line  segment  formed  by  the  intersection  of 
the  oath  with  the  set  of  SCENARIO  grid  cells  through  which  it  passes. 
After  initializing  each  integral  register  to  zero,  RFPROP  begins  to  loop 
through  the  set  of  data  points.  At  each  point,  the  submodule  DATAPT 
(described  in  Section  3.4)  is  called  to  retrieve  the  set  of  environmental 
data  for  the  eight  grid  points  surrounding  the  point  and  to  interpolate 
these  to  the  point.  The  environmental  data  set  consists  of  the  following: 
(1)  the  mean  electron  density  [cm""],  (2)  the  electron  temperature  [°K], 
(3)  the  electron  density  standard  deviation  about  the  mean  [cm"3],  (4)  the 
rms  electron  density  [cm-*],  (5)  the  electron-neutral  collision  frequency 
[sec-1],  (6)  the  set  of  electron-ion  collision  freguencies  [sec"1]  (the 

electron-ion  collision  freguency  is  a  function  of  carrier  frequency),  (7) 
the  three  components  of  the  ion  velocity  relative  to  the  tangent  plane  of 

the  point  [cm  sec'1],  (8)  the  partial  time  derivative  of  the  mean  electron 

density  [cm" "sec" 1 ],  (9)  the  outer  striation  scale  sizes  in  two  directions 
perpendicular  to  the  field  [cm]  (currently,  these  two  values  are  equal), 
(10)  the  orientation  angle  of  the  major  transverse  striation  dimension 

(currently  zero),  (11)  the  outer  striation  scale  size  parallel  to  the 
field  [cm],  (12)  the  isotropic  inner  striation  scale  size  [cm],  (13)  the 

three  components  of  the  dipole  magnetic  field  relative  to  the  tangent 

piano  of  the  point  [Gauss],  (14)  the  three  components  of  the  mean  electron 
density  qradient  relative  to  the  tangent  plane  [cm-4],  and  (15)  the  six 

components  of  the  second  spatial  derivative  of  the  mean  electron  density 
relative  to  the  tangent  plane  [cm-^].  The  components  of  the  vectors  and 
the  second  derivative  are  then  transformed  to  a  path  oriented  coordinate 

system  as  described  later. 
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2.2.1 


Geometry  of  Propagation  Path. 


All  calculations  are  made  in  (or  transformed  to)  the  coordinate 
system  aligned  with  the  line  of  sight  (LOS)  along  the  path.  The  orienta¬ 
tions  of  the  coordinate  axes  of  this  system  are  as  follows:  the  W-axis  is 
oriented  along  the  line  of  sight  pointing  away  from  the  transmitter  (the 
path  origin)  toward  the  receiver,  the  U-axis  is  perpend  icu  1  ar  to  the  path 
and  lies  in  the  "vertical"  plane  containing  the  path  and  the  center  of  the 
earth,  and  the  V-axis  is  perpend  i  cu  1  ar  to  the  path  and  perpendicu  lar  to 
the  "vertical"  plane.  Note  that  the  U-axis  always  points  away  from  the 
center  of  the  earth.  The  transformation  from  tangent  plane  coordinates  to 
LOS  coordinates  is  defined  by  a  clockwise  rotation  of  (tt/2+x)  about  the  Z 
(vertical)  axis  followed  by  a  clockwise  rotation  of  (n/2-e)  about  the 
Y-axis  where  is  the  elevation  of  the  path  measured  from  the  transmitter 
foriqin)  to  the  receiver,  and  X  is  the  azimuth  of  the  path. 

In  addition  to  the  tangent  plane  (X-Y-Z)  and  LOS  (U-V-W)  coordi¬ 
nate  systems,  two  other  systems  are  used  in  RFPROP  (see  figure  2).  These 
are  the  field  al igned  coordinate  system  (r-s-t)  in  which  the  power  spec¬ 
tral  density  is  expressed  in  Equation  2-1  and  a  system  intermediate 
between  the  field  aligned  (r-s-t)  system  and  the  LOS  (U-V-W)  system  which 
Wittwer  designates  as  the  (x-y-z)  system.  Note  that  this  system  must  NOT 
be  confused  with  the  tangent  plane  (X-Y-Z)  system.  In  the  field  aligned 
system  the  t-axis  is  parallel  to  the  field  vector  and  the  r-axis  is 
oriented  in  the  plane  perpend icul ar  to  the  field  vector  along  the  direc¬ 
tion  in  which  the  scale  size  is  maximum.  (The  s-axis  points  in  the  direc¬ 
tion  along  which  the  scale  size  is  minimum.)  Currently,  the  SCENARIO 
model  provides  only  one  scale  size  perpend icu 1 ar  to  the  field  direction, 
and  the  choice  of  the  r  and  s-axes  is  arbitrary.  Until  the  orientation  of 
the  r-axis  is  specified  explicitly  by  SCENARIO,  we  shall  assume  that  the 
r-axis  and  the  s-axis  are  equivalent.  The  x-axis  of  the  (x-y-z)  system  is 
perpend  icu  1  ar  to  both  the  path  and  the  field  direction,  that  is 
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Figure  2.  Geometric  transformation  from  field  oriented  coordinates  to  line 
of  sight  (LOS)  coordinates.  In  (a)  the  magnetic  field  is 
oriented  along  t,  and  the  LOS  lies  along  z.  Rotations  about  two 
axes  suffice  to  align  t  and  z.  In  (b)  R  points  toward  the 
center  of  the  earth,  and  together  with  z,  it  defines  the 
■vertical"  plane  containing  the  LOS.  The  orientation  of  this 
plane  is  invariant  along  the  LOS,  and  the  components  of  all 
vectors  transverse  to  the  LOS  are  specified  relative  to  it. 
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X 


t  X  z 
sin  $ 


(2-4) 


where  4>  is  the  angle  between  the  path  and  the  field  vector 

cos$  =  t  •  z  =  B  / | B |  (2-5) 

and  Bz  is  the  field  component  along  the  path.  More  briefly,  the  y-axis  is 
oriented  along  the  field  component  in  the  plane  perpendicular  to  the  path. 

The  transformation  from  (r-s-t)  coordinates  to  (x-y-z)  is  defin¬ 
ed  by  two  rotations:  first,  the  (r-s-t)  system  is  rotated  counterclock¬ 
wise  about  the  t-axis  through  an  angle  to  align  the  r-axis  with  the 
x-axis;  next,  the  system  is  rotated  counterclockwise  about  the  x-axis 
through  the  angle  T  to  align  the  t-axis  with  the  z-axis.  Note  that  if  the 
plasma  structure  is  isotropic  in  the  r-s  plane,  we  may  set  <t>  =  0.  One 

additional  rotation  is  required  to  align  the  (x-y-z)  system  with  the 

(U-V-W)  system.  To  align  the  x-axis  with  the  U-axis  (in  the  vertical 

plane),  rotate  the  (x-y-z)  system  counterclockwise  about  the  z  =  W-axis 
through  the  angle  0,  where 

tanO  =  Bu/Bv  (2-6) 

Note  that  the  W  and  z-axes  coincide.  Here,  we  often  refer  to  Wittwer's 

z-axis  as  the  W-axis.  In  contrast  to  the  three  other  systems,  the 
orientation  of  the  (U-V-W)  is  fixed  in  space  along  the  entire  length  of 
the  path.  All  integrals  are  calculated  in  this  fixed  system. 
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Table  2.  Summary  of  coordinate  transformations  for  the  four  systems  used 
toy  RFPROP  module. 


T  anqent 

P 1  ane 

to  LOS.  . 

(X-Y-Z) 

to 

(U-V-W) 

(1) 

clockwise  ir/2  + 

X 

about 

Z 

(2) 

clockwise  n/2  - 

£ 

about 

Y 

Magnetic 

to  Intermediate 

LOS 

(r-s-t) 

to 

(x-y-z) 

(1) 

counterc lockwise 

about 

t 

(2) 

counterc lockwi se 

$ 

about 

X 

Intermed i ate 

LOS  to  LOS. 

(x-y-z) 

to 

(U-V-W) 

(1) 

counterclockwise 

9 

about 

W  =  z 

To  simplify  evaluation  of  the  PSD  and  the  signal  parameters, 
dpfine  the  "scale  size  operator"  in  magnetic  coordinates  as  the  diagonal 
matrix: 

L  =  d iaq(L 2,  L  2,  L2)  [cm2]  (2-7) 

To  express  the  PSD  function  in  the  LOS  coordinates,  we  transform 
the  elements  of  the  scale  size  operator  by  the  second  pair  of  rotations  in 
Table  2 : 

L2  =  L'V.os2*  +  [2sin2^  Frm2!  f  ?-Rl 

XT  S  u  -*  k  ' 


L2  =  (L2sin2-t>  +  L  2cos  7j>)cos  ?  +  L2sin2f  [cm2]  (2-9) 
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(2-10) 


L?  =  (L?sin?+.  +  L?cos  ?f)s  in  +  L,?cos?* 


L  =  (L‘ 
xy 


Lr)cos ^cos^s  in 


L  =  ( L ' 

V  ~T  '  W 


'xz 


L^)s  intcosfs  in* 


L  =  (L  2  -  L'sin2*  -  L  2cos  2  t)cos  ts  in  $ 
yz  t  r  s 


[cm2] 

[cm ?]  (2-11) 

[cm2]  (2-12) 

[cm2]  (2-13) 


(See  Equation  5-a  through  5-f  of  DNA-IR-82-022 . )  And,  finally,  to  express 
the  transverse  elements  of  the  operator  in  the  fixed  (U-V-W)  system: 


L2  -- 

L  2cos 

o  7.7 

■n  +  L  Sinn  - 

2L  coses  inn 

[cm  ?] 

(2-14) 

u 

y 

X 

xy 

LV?  " 

L2cos 

+  LySin?o  + 

2Lv..cos  esin  n 

[cm2] 

(2-15) 

Luv 

=  (L2  . 

y 

-  L2)cosqsinn  - 

•  L  (  sin ?n-cos 2n) 
xyv  ' 

[cm2] 

(2-16) 

(See  equations  preceding  Equation  18  of  DNA-IR-82-022.)  Note  that  Wittwer 
has  switched  the  sense  of  the  U  and  V-subscripts  of  the  scale  size  opera¬ 
tor  in  order  to  simplify  the  form  of  some  of  the  equations  which  follow. 


The  "effective"  structure  scale  size  along  the  path  (W-axis)  is 
defined  to  be 


'/7 


LrLsLt 

v'det(L) 


where 


det<L>  =  LuLv  -  <Luv)?  =  LxLy  - 


[cm]  (2-17) 


[cm4]  (2-18) 
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is  the  determinant  of  the  scale  size  operator  on  the  U-V  plane  --  or  on 
the  x-y  plane.  (See  Equation  10  of  DNA-IR-82-022 . ) 

The  total  time  derivative  of  the  mean  electron  density  is  com¬ 
puted  as  the  sum  of  the  partial  time  derivative  (as  computed  by  SCENARIO) 
and  the  scalar  product  of  the  density  gradient  and  the  velocity  of  the 
point  on  the  path 

dN  9N 

— -  =  — -  +  V.  *qrad(N  )  [cm_3sec_1]  (2-19) 

dt  at  Pl  e 


where 


V  + 
org 


zend 


-V  ) 
org7 


[cm  sec-1]  (2-20) 


and  V  and  V  ,  are  the  velocities  of  the  transmitter  and  receiver, 
orq  end 

respect  i  vely,  and  z  and  z0nd  are  the  distances  to  the  point  and  to  the 
receiver  from  the  transmitter,  respectively.  The  second  and  third  time 
derivatives  are  estimated  by  assuming  a  power  law  function  for  the  mean 
electron  density: 


Ne(t )  =  at 


[cm-3]  (2-21) 


Then 


and 


d2N 


dt' 


1 


dNQ  - 

-  f-^)2 

Ne  dt 


1 


dN. 


(— ) 
t  dt 


[cm- ^sec" 2] 


(2-22) 


d3N 


dt ' 


-  i-  f 


d2N, 


d2N  0  ,  dNQ 

— -f  -  -  (— )  ( — 2) 

N0  dt?  N0  dt 


dt' 


[cm" 3sec- 3] 


(2-23) 
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Next,  the  components  of  the  weighted  effective  plasma  velocity 
at  the  point  are  calculated 

*eff  '  f— "-1  f*Vv>  1“  (2-24) 

where  V-Qn  is  the  local  plasma  velocity  relative  to  a  fixed  coordinate 
system,  and  Vpt  is  the  velocity  of  the  point  (Equation  2-20).  (See  equa¬ 
tions  following  Equation  14  of  DNA-IR-82-022 .)  This  is  used  to  compute 
the  quantity  which  in  turn  is  used  to  compute  the  perpendicular  signal 
decorre 1  at  ion  time: 


2  _ 
Dv  - 


LV 
y  x 


-  2L  v  v 
xy  x  y 


+  lV 

x  y 


det(L) 


[sec-2]  (2-25) 


(See  equation  followinq  Equation  14  of  DNA-IR-82-022.)  Although  we  have 
chosen  to  calculate  it  in  the  (x-y-z)  system,  it  is  a  scalar  quantity  and 
may  be  calculated  in  any  coordinate  system. 

Next,  calculate  Wittwer's  Bp  coefficient  --  the  coefficient  of 
the  first  term  in  the  Taylor  series  expansion  of  the  differential  decorre¬ 
lation  function.  (See  Equation  B-l  of  DNA-IR-82-022.) 


,rB  n(o2)W? 


(0?te2)(n-1)/?  ) 
Cn2n-2r(n-l) 


(2-26) 


where  n  is  a  unitless  quantity 


o 


? 


2,V 


2L  xy  +  L2x2 
xy  J  y 

det(L) 


(2-27) 


and  r  is  the  ratio  of  the  inner  to  the  (minimum)  outer  scale  size,  with 
Eouation  2-26  valid  in  the  range 
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(2-28) 


m  minimum  (2,2n-2) 


(2-29) 


.-_.Qii.Li 

?n'?;'(n-i) 


(2-30) 


where  K  ,  is  the  modified  Bessel  function  of  order  n-1.  For  the 

n- 1 

range  of  values  of  .  considered  here,  it  is  easy  to  show  that  c^  does 
not  differ  significantly  from  1.  (See  Fquation  9.6.9  and  Table  9.8  of 
Abramowit/  and  Stegun  '. ) 


The  coefficient  R^  is  then  approximated  as 


R  ■  —  —  f(n)h  (t) 

n  ('n  -  ;•( n-1) 


(2-31) 


(n-1) 

(n-1) 


-  ; .  m 


1.2  ■  n  v  2.0 


2.0  ■  n  •  4.0 


(2-32) 


\,<  > 


4-?n 
1/3  -  I  n( ■  ) 


1/3  ♦ 


1.2  >  n  2.0 


n  2.0 


2.0  •  n  -  4.0 


(2-33) 
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(See  Appendix  B  of  DNA- 1 R-82-02 2 . )  Note  that  both  Bn  and  p  are  unit¬ 

less  and  independent  of  the  propagation  frequency. 

2.2.2  Point-Wise  Integration  Along  Propagation  Path. 

RFPROP  now  begins  to  loop  over  the  set  of  simulation  frequen¬ 
cies,  to  calculate  each  integrand  and  sum  that  integrand  into  the  appro¬ 
priate  primitive  integral  reqister.  For  these  calculations,  we  need  the 
frequency  related  parameters 


angular  frequency: 

u)  = 

2,fc 

[rad  sec-1] 

(2-34) 

wave  length : 

x  = 

c/fc 

[cm] 

(2-35) 

wave-number : 

K  = 

2  tt/  X 

[rad  cm-1] 

(2-36) 

and  the  unitless  variance  of  the  index  of  refraction  of  the  signal 


(2-37) 


where  f  is  the  signal  carrier  frequency  (Hz),  c  is  the  speed  of  light 
(2.99RE+10  cm  sec'1),  re  is  the  classical  radius  of  the  electron 

p 

(2.8184E-13  cm),  and  is  the  variance  of  the  density  at  the  point  as 
calculated  by  SCENARIO.  In  order  to  eliminate  unnecessary  computuat ions , 
all  calculations  at  a  point  are  omitted  if  both  the  local  mean  electron 
density  and  its  standard  deviation  fall  below  an  arbitrary  threshold  value 
(currently  set  at  1  cm"3  for  both  parameters). 
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The  differential  phase  variance  at  the  point  is 


? 


w  - Y  k2  "1 


[rad 2cm_  !] 


(2-38) 


where  is  the  effective  scale  size  along  the  LOS  (Equation  2-17),  and  y 
is  a  (unitless)  constant 


2/?  r(n-l) 
Y  "  r(n-3/2) 


for 


n  >  1.5 


(2-39) 


(See  Equations  8  and  9  of  ONA- IR-82-02 2 . )  The  total  phase  variance  along 
the  path  is  obtained  by  multiplying  the  differential  phase  variance  by  the 
path  segment  length  at  the  point  and  summing  over  the  N  points  along  the 
path 


a 


2 

4> 


N 

y 

i 


[rad2]  (2-40) 


Also,  calculate  the  phase  weighted  sum  of  the  Bn  coefficient 
(Equation  2-31)  which  is  used  to  generalize  the  frequency  selective  band¬ 
width  calculation  from  a  delta  layer  approximation  to  a  layer  of  finite 
thickness 

N  ,2 

sum(h)  =  V  B  —  az  [tad2]  (2-41) 

i  n  dz 


(See  equation  preceding  Equation  32  in  DNA-IR-82-022 .) 

Next,  we  calculate  the  [unitless]  function 
.  <LxtL(>(zend-z>(z/zend> 

M(z)  2 — rein - 


(2-42) 


and  the  three  associated  primitive  integrals  required  to  compute  the  mean 
square  log  amplititude  fluctuation  of  the  signal 
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(2-43) 


where 


f(n,z) 


and 


where 


N  o 

sum(y0)  =  —  Y  az  (n-l)  f(n*z)  ^q(z) 

2.  \  dz 

sum(Xl)  =  1  V  az  (n-1)  f(n,z)  I^z)  (2-44) 

sum(x2)  =  —  y  AZ  (n-1)  f(n,z)  ^(z)  (2-45) 


f'(n)  +  -  M(z)  > 

=  ( _ - _ )  {1  -  —  exp  [-  (^il 

1  1  +  M( z)  j  12  3 


10M(z) 


-  I)2]! 


f 1 ( n)  =  1.1  -  maximum  (0,  (n-2.4)/2) 


(2-46) 

(2-47) 


In(z)  = 


2n,  6-2n  ,6-2n^  M(z)n"1 


i6crv2 


16cf n  In  (c2/ai)  M(z) 


6-2n 

n-1 


n*3 

n=3 


Ix(z)  =  8[I  aj  +  cj  ( a2~c2) ]  M(z)-1 


2  Z'  "  3C1  Sa2 


(2-48) 


(2-49) 

(2-50) 


al 

=  minimum  fc2,Cj/M(z)] 

(2-51) 

a2 

=  maximum  [ c2 ,c1/M( z) ] 

(2-52) 
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with  cx  =  0.5  and  c2  =  0.84.  (See  equations  following  Equation  A-3  of 
DNA-IR-32-022 . )  Note  that  M  and  all  the  primitive  sums  are  unitless  quan¬ 
tities.  This  derivation  for  the  mean  square  log  amplitude  is  valid  for 
all  values  of  M  and  for  n  within  the  range  1.2  _<  n  _<  4.  Also,  RFPROP 
stores  the  [unitless]  auxiliary  sums 


i  N  ,2 
1  c»  4  _  d  G 


2  i  dz 


az  [M (  z )  ]  2 


(2-53) 


I  V  Az  M(  z)  n_1 


(2-54) 


-  -  AZ 


ln[M(z) 1M( z) 2  n=3 


which  will  be  used  to  evaluate  the  Rayleigh  phase  variance, 
tions  A-8a  through  A-8c  of  DNA-IR-82-022 . ) 


(See  Equa- 


Calculate  the  integral  to  be  used  to  compute  the  signal 
decorrel ation  time  in  the  direction  perpendicul ar  to  the  LOS 


sum(r)  =  ;  AZ  JL-  Bn  (ovz/zend)r 


[ sec~m] 


(2-55) 


where  8n  is  the  coefficient  from  Equation  2-32  and  ov  is  given  by  Equation 


2-25.  (See  Equation  14  of  DNA-IR-82-022.) 


Sum  the  weiqhted  plasma  velocity  relative  to  the  receiver  in  the 
directions  parallel  and  perpend icul ar  to  the  LOS: 


Sl,n(v;)  •  '  >.z  i£  l(Vion-Vend)-  1 1  [cm  sec-’-] 


(2-56) 


sumfvjJ  =  v  Az  -  V 


trans 


[cm  sec' (2-57) 


where 


trans 


i r rv -  -v  ,  Vu l2  +  r ( v .  -v  .vl2}1/2 

ion  end  !  '  ion  end ' 


[cm  sec'1]  (2-58) 


Calculate  and  sum  the  increments  of  the  primitive  integrals 
needed  to  calculate  the  moments  of  the  coherence  function  at  the  receiver: 


sum(Cu) 


=  V  B  az  r(z/zp  ,)2  lm/2 

j  n  dz  enc*  det(L) 


[cm_m]  (2-59) 


sum(C„)  =  Y  B„  —  az  f (z/z  .) 


2  Lv  im/2 


end'  det(L) 


[cm'm]  (2-60) 


SM(cu»>  '  [  Bn  ST  Az  ,(z/2end»2  Siry,,n/2  s'5"ad»)  [cm"m] 

(See  Equations  20-a  through  20-d  of  DNA-1R-82-02 2. )  The 
quantities  are  computed  at  the  transmitter  by  replacing  z  in 
2-59  through  2-61  by  (zen(j  -  z). 

The  signal  time  delay  variance  (o2)  and  frequency 
bandwidth  (f0)  have  been  derived  by  Wittwer  (see  Equations 
DNA- IR-82-022 ) 


(2-61) 

analogous 

Equations 

selective 
31-33  of 


_  0  H  zend  .  z  .  ,  0 

c  K  K2  o  z2  o  z’2  u 

+  I2(z‘)  +  21 2  (z)  ]) 


[sec2]  (2-62) 


where 


is  the  phase  variance  of  the  Rayleigh  distributed  com¬ 
ponent  of  the  signal 

is  the  parameter  derived  in  Equation  2-127 
z end  1S  the  length  of  the  path  [cm] 

and 

2 '  da^ 

1,(2')  =  2/  B  dz  (_*]  z2L2|det(L)  I" 1  (2-64) 

u  0  n  dz  u 

z'  da2 

I  (z1)  =  2/  Bndz  (-_*)  z 2L 2 f det ( L )  1"  1  (2-65) 

o  dz 

z‘  da2 

W2'*  =  Bndz  z  Luv(det(L)  1‘  [2-bb) 

o  az 

where  (da2/dz)  is  the  differential  phase  variance  (Equation  2-38),  B  is 

<p  II 

defined  in  Equation  2-31,  det(L)  is  the  determinant  of  the  scale  size 
operator  (Equation  2-18),  and  L2  etc.,  are  the  transformed  elements  of  the 
scale  size  operator  (Equations  2-14  through  2-16). 

Since  the  calculation  depends  on  a  linear  combination  of  the 
sauares  of  the  innermost  integrals  (Equations  2-64  through  2-66),  we  shall 
limit  the  remainder  of  this  discussion  to  the  evaluation  of  the  single 
3-level  nested  integral 
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(2-67) 


r  *  Xrend  xf  iL  I/dc  e3  f(c)l2  [cm-3] 
0  x2  x' 2  0 


All  the  other  integrals  computed  by  RFPROP  are  evaluated  using  a  mean- 
rectangular  integration  algorithm  (see  Equation  2-3).  The  accuracy  of 
this  technique  is  unsatisfactory  for  the  triply  nested  integral  of  Equa¬ 
tion  2-67,  however.  First,  we  simplify  the  integral  using  integration  by 
parts  of  the  outermost  integral  to  yield  the  2-level  nested  integral 


=  fend  ^  fl  -  _L_)  ff  dx'x,2f(x')  l2  [cm*2] 
o  x2  ^  xend  o 


(2-68) 


Then,  approximate  the  integrand  f(x)  by  a  piecewise  linear  function 
between  the  n+1  discrete  points  on  the  path 


04X 

+ 

Bl 

on 

x  0  < 

X 

< 

*1 

a2x 

+ 

62 

on 

Xj  < 

X 

< 

*2 

anx 

+ 

6n 

on 

xn-l 

< 

X 

—  xn 

(2-69) 


where 


"j  =  (yj-yj-i)/(xj-xj-i} 


(2-70) 


"j  =  (xjyj-r 

f or  j  =  1 ,  2 ,  .  .  . ,  n  . 


(2-71) 


It  can  be  shown  that  the  innermost  integral  can  be  calculated  by 
the  recursion  relation 
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where 


Also, 


Then, 


where 
r  ( i ) 

j 


j  =  eRi 4  r(3>uj)ay(j)  - 1  r(2)(V  <2-72> 


6  .  r  x  .  ,  /x  . 
J  J-l  J 


r  ^ n  ^  ( F,  . )  =  1  +  F,  .  +  £ 2  + 
y  j  ,i 

Ay( 0  =  y .  -  y .  . 

J  J  J-l 


-  + 


Ayvy  =  c  ,-ys  -  y. 


■J'j  yj-l 


define  the  quantity 


Cl  --  c3  [yl  ,  -  ( 1-e . ) ' L  (1  C  .  Ay(  \ ) 
J  J  1  J-l  J  M  J  J 


(2-73) 

(2-74) 

(2-75) 

(2-76) 


I&y(2))]  (2-77) 


it  can  he  shown  that  the  integral  of  Equation  2-68  rnay  be  written 


=  r( 1 )  - 


end 


r( 1  ■ 

n 


[cm'2] 


(2-78) 


the  r.  are  calculated  by  the  recursion  relations 
J 


:  r 


l\h  +  X4  I  (  1-r  .  )-  1  fl_  r(5)(r  .)(ay(l))2  +  1_  r(3)(r  ,)(av(2.))2 
J-l  J  '  1  "96  j  3  J  36  V  V  '  J  1 


"96 

i—  r  ( 4  )( r,  . )  Ay(  0  Ay  (2. )  1 

30  J  J  i  ' 


+  Cl  [i  r  (  1  )  (  r  . )  Ay(l)  -  |  Ay(2.)  -  1  Cl  (l-S-'2)]}  (2-79) 
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and 


r(  ?)  =  r(  ?),  +  xs  ll-r . )_  1  —  r(  6)  ( 5 . )  (  Ay(  .0  )2  +  i_  r( 4)  (?,.)(  av(  ?) )  2 

J  J-l  J  J  11?  J  J  45  J  '  J 

-  L  r(5)(r.  )&y(!)Ay(2)l 

36  j  j  j 

+  C^l  r(2)(5.)Ay(|)  -  I  r(0(C  )Ay(?)  -CUl-q1)]}  (2-80) 

J  6  J  J  3  JJ  J  J  1 

The  summation  registers  containing  the  sums  of  the  terms  of 
Equations  2-79  and  2-80  are  incremented  inside  the  grid  only.  When  the 
LOS  extends  beyond  the  edge  of  the  grid  (i.e.,  whenever  the  receiver  is 
not  located  inside  the  grid),  it  is  necessary  to  correct  the  sum  to  ac¬ 
count  for  the  finite  distance  from  the  edge  of  the  grid  to  the  receiver, 
rt  can  be  shown  that  the  necessary  correction  factor  to  be  added  to  Equa¬ 
tion  2-78  is 


e 


-  1)?  (  "  )4 
end 


(2-81) 


where  xn  is  the  distance  from  the  origin  ( transmi tter )  to  the  far  edge 
of  the  grid,  x  ^  is  the  length  of  the  LOS,  and  is  the  innermost  inte¬ 
gral  obtained  from  Equation  2-72.  (For  a  more  complete  explanation  of  the 
derivation  of  these  equations,  refer  to  MRC-N-273 16 . ) 


All  the  remaining  basic  parameters  calculated  by  RFPROP  are  used 
to  derive  the  mean  ionization  effects.  We  calculate  the  real  part  of  the 
index  of  refraction,  the  absorption  coefficient,  and  the  effective 
electron  density  using  the  App 1 eton-Hartree  equations.  PRPSIM  considers 
only  the  contribution  of  electron-ion  collisions  and  el ectron-neutral 
collisions.  The  magnetic  field  and  ion-neutral  collisions  can  be 
neglected  at  the  propagation  frequencies  of  interest  here.  Now  let 
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X 


(2-82) 


enm 


3 . 1833x10 9  Ne/u>2 


where  e  and  m  are  the  electronic  charge  (4.803E-10  statcou lombs)  and  mass 
(9.107E-28  gm) ,  respectively,  is  the  permittivity  of  free  space 
f(4rr)_  1  =  7.958E-2  in  cgs  units),  is  the  local  electron  density 
[cm-3],  and  u>  is  the  angular  frequency  of  the  signal  [rad  sec-1];  and  let 

v  .  +  V 

Z  =  - —  (2-83) 

/Jk) 


where  v  .  and  v  are  the  electron-ion  and  electron-neutral  collision  fre- 
c  1  en 

quencies,  respectively.  These  parameters  are  then  modified  to  account  for 
the  velocity  dependence  of  the  electron  collision  frequency: 


X 


eff 


(2-84) 


where 


-eff 


G  =  i 


H 

l+Z/2.5 

l+Z/1.5 

1 


V  >  V  . 


en  ei 


<  V  • 

en  —  ei 


(2-85) 


(2-86) 


and 


H  H 


1  +  .  15Z 
1  +  .05Z 

1 


>  v 
en  e 


v 


ei 


(2-87) 


(See  Shkarofsky,  1961 4 . )  The  real  part  of  the  [unitless]  index  of  refrac¬ 
tion  is  given  by 
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where 


1 


=  -W(l -Y)  +  /(1-Y)2  +  Z2  .  Y2 


n 


eff 


Y  = 


eff 


1  +  Z 


eff 


(See  Equation  2.93  of  Davies,  1965s.) 


(2-88) 


(2-89) 


The  linear  coefficient  of  absorption  is  given  by 


20  1  U  y? 

’  ln(10)  2c  “  eff 


1.449x10“  10  J? 

u 


YZ 


eff 


[dB  cm"  l] 


(2-90) 


(See  Equation  2.98  of  Davies,  1965s.)  If  one  can  assume  that  the  carrier 
frequency  is  large  compared  to  the  collision  frequencies,  the  coefficient 
of  absorption  may  be  approximated  as 

a  =  Ne(vei+  ven)  [dBcm-1]  (2-91) 

yu) 

RFPROP  calculates  the  mean  absorption  coefficient  for  the  cell 

00 

a  =  /  cc(N  )  p(N  )  dN  (2-92) 

o  K 

where  p(N  )  is  the  unknown  distribution  function  of  the  electron  density 
within  the  cell.  From  Equations  3-28  and  3-30,  it  is  seen  that  the 
neutral-electron  collision  frequency  (per  electron)  is  independent  of  the 
electron  density,  and  that  the  ion-electron  collision  freouency  is  propor¬ 
tional  to  the  electron  density.  Thus 
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a  =  /  (a  jNj  +  a?Ne)  p(N  )  dN0 

o 


aiNe  +  a2Ne 


(2-93) 


Hence,  the  mean  absorption  coefficient  can  be  computed  as  the  sum  of  two 
coefficients  with  the  first  being  a  function  of  the  ion-electron  collision 
frequency  and  the  root-mean-square  electron  density,  ana  the  second  being 
a  function  of  the  neutral -e 1 ectron  collision  frequency  and  the  mean  elec¬ 
tron  density. 


The  effective  mean  electron  density  at  the  signal  propagation 
frequency  is  given  by 


Ne 


eff 


1  +  Z 


eff 


[cm*3] 


(2-94) 


where  Ng  is  the  mean  electron  density  provided  by  DATAPT,  and  Zg^  is  the 
Z  coefficient  calculated  using  only  the  electron-neutral  collisions.  Note 
that  the  effective  density  differs  significantly  from  the  actual  density 
only  when  the  signal  propagation  frequency  is  comparable  to  or  smaller 
than  the  electron-neutral  collision  frequency.  This  value  is  added  to  the 
electron  content  register  to  calculate  the  mean  total  electron  content 
(TFC)  alonq  the  path. 

N 

TEC ■  =  >  A2  Nepff  [cm'2]  (2-95) 


The  first  three  time  derivatives  of  the  local  mean  electron 
density  (calculated  by  SCENARIO  and  in  Equations  2-22  and  2-23  above, 
respectively)  are  also  summed  into  individual  registers  to  calculate  the 
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three  time  derivatives  of  the  "TEC>.  Similarly,  the  total  absorption  is 
calculated  as  the  sum  of  the  absorption  coefficients  (Equation  2-93) 
multiplied  by  the  path  segment  lengths. 


The  incremental  Faraday  rotation  of  the  plane  of  polarization  of 

the  electric  vector  of  the  signal  is  calculated  for  the  case  of  quasi¬ 

longitudinal  propagation  i.i  the  absence  of  absorption: 

2  MB  N  B 

—  =  <  — _ U-T  J e L— -  =  2.365*104  Jl-DL  [rad  cm’1]  (?-96) 

dz  4  n 2  r  0  m  m  f2  f2 

where  c  is  the  speed  of  light,  e0  is  the  permittivity  of  free  space,  e  and 
m  are  the  charqe  and  mass  of  the  electron,  respectively,  Ng  is  the 
effective  mean  electron  density  [cm-3],  and  B w  is  the  component  of  the 
magnetic  field  [Gauss]  parallel  to  the  LOS.  (See  Equation  4.39  of  Davies, 
1965  s . )  The  total  angle  of  rotation  is  calculated  as  the  sum  of  the 
incremental  rotations  multiplied  by  the  path  segment  lengths. 

The  last  frequency  dependent  inteqrands  to  be  computed  are  the 
weighted  distance  to  the  point  and  the  weighted  square  of  the  distance 
with  the  phase  variance  as  the  weighting  function.  These  are  used  to 
provide  measures  of  the  distance  to  the  most  intense  region  along  the  path 
and  of  the  extent  of  the  active  region.  Let 

N  2 

sum(Z)  =  "  z  az  (——)  [cm]  (2-97) 

i  dz 

’1  .  ? 

sumfZ7)  -  "  z2  tz  r_L)  [cm2]  (2-98) 

i  dz 

This  point,  marks  the  end  of  the  loop  for  which  the  inteqrands 
arc  calculated  as  functions  of  frequency  and  also  the  end  of  the  loop  over 
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the  points  along  the  path  at  which  the  integrands  are  summed  into  their 
appropriate  integral  registers. 


2.2.3  Calculation  of  Signal  Parameters  from  Primitive  Integrals. 


Next,  RFPROP  begins  to  loop  over  frequency  to  derive  the  signal 
parameters  from  the  primitive  integrals  which  were  summed  over  the 
individual  points  along  the  path.  First,  the  principal  moments  and  the 
orientation  of  the  decorrelation  matrix  at  the  receiver  are  computed: 


=  r  sum(Cu)  ] 2/m 

[cm-2] 

(2-99) 

=  f  sum(Cv)  I2/"1 

[cm-2] 

(2-100) 

=  sign[ sum(Cuv)]|sum(Cuv) |2/m 

[cm-2] 

(2-101) 

where  sum(Cu),  sum(Cv),  and  sum(Cuv  )  are  the  primitive  integrals  by  Equa¬ 
tions  2-59  through  2-61.  (See  Equations  20a-d  of  DNA-IR-82-022.)  Find 

the  maximum  and  mininum  eigenvalues  of  the  matrix,  C  and  C  : 

P  q 


=  1  ffC  +C  ) 

+  /(C  -C  )2  +  4C2  1 
u  V '  uv  1 

[cm-2 ] 

(2-102) 

=  I  [fc  +C  ) 
2  u  v ' 

-  /fc  -C  )2  +  4C2  ] 

1  U  V  UVJ 

[cm-2] 

(2-103) 

and  the  orientation  of  the  maximum  eigenvector: 


£ 


A  arctan 
2 


2C 

[— ! ^-] 

C  -C 
u  v 


[rad]  (2-104) 
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(See  equations  preceding  Equation  21  of  DNA-IR-82-02 2.)  Note  that  Cp 
and  denote  the  maximum  and  minimum  eigenvalues,  respectively  of  the 
decorrel at  ion  matrix.  The  analogous  decorrelation  matrix  at  the  transmit¬ 
ter  is  calculated  also.  The  decorrelation  distances  along  the  p  and  q 
axes  are  derived  from  the  eigenvalues  of  the  decorrel at  ion  matrix: 

B(n,o2) 

i  =  _ _l_  [cm]  (2-105) 

p  /C 

P 

B(n,a2) 

l  =  _ _JL  [cm]  (2-106) 

P  Jr 


where  B(n,o2)  is  the  empirical  function  defined  in  Equation  2-112.  (See 
Equation  41  of  DNA-5304D1.)  Note  that  l  and  i  denote  the  minimum  and 

r  M 

maximum  decorrelation  distances,  respectively.  The  decorrelation  distance 
is  defined  as  the  distance  from  the  LOS  on  the  perpendicular  plane  at 
which  the  mutual  coherence  function  decreases  by  a  factor  of  2.718.  (See 
Equation  21  of  DNA-IR-82-022.)  The  mean  square  angles  of  the  signal 
energy  arrival  at  the  receiver  and  transmitter  are  derived  from  the  decor¬ 
relation  matrix  also: 


2Cp 

P  K2 


[  r  ad 2  ] 


(2-107) 


c 


? 

q 


K2 


[rad2] 


(2-108) 
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(See  Equations  26a  and  26b  of  DNA-IR-82-022.)  Note  that  a2  and  a2  denote 
the  maximum  and  minimum  values,  respectively.  Also  note  that  each  of  the 
Equations  2-99  through  2-108  denotes  the  calculation  of  two  quantities  -- 
one  at  the  receiver  and  one  at  the  transmitter. 


The  signal  decorrelation  time  at  the  receiver  in  the  direction 
parallel  to  the  LOS  is  approximated  by 


T  = 


3.5  K 


r c p / 3  +  c2/3!3/2  <v„> 


[sec] 


(2-109) 


where  <v  >  is  the  weighted  average  magnitude  of  the  plasma  velocity 
relative  to  the  receiver  and  parallel  to  the  LOS  (from  Equation  2-56) 

<vn>  =  sumfv^/a2  [cm  sec-1]  (2-110) 


(See  Equation  2-56  above  and  Equation  22  of  DNA-IR-82-022.) 


The  signal  decorrelation  time  at  the  receiver  in  the  plane 
perpendicular  to  the  LOS  is 


T1  = 


rsum( x) 


[sec] 


(2-111) 


where  sum( t)  is  the  primitive  integral  of  Equation  2-55,  m  is  defined  in 
Equation  2-29,  and  B  is  a  [unitless]  function  of  the  spectral  index  and 
the  total  phase  variance: 

B(n,a2)  =  min [ 1 , (-0. 34n  2  +  2.51n  -  2.0) ] fsum(h)  1 1 /m  (2-112) 


where  sum(h)  is  the  phase  weighted  sum  (Equation  2-41)  of  the  Bn 
coefficient  calculated  in  Equation  2-31.  (See  Equation  14  of  DNA-IR-82- 
022.)  The  B  function  imposes  an  emperically  determined  lower  limit  for 
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the  perpend i cul ar  decorrel at  ion  time  when  the  total  phase  variance  along 
the  path  becomes  vanishingly  small. 

The  smaller  value  of  the  tJ_  -  t  pair  is  designated  as  the  sig¬ 
nal  decorrelation  time: 

t 0  =  min(xj_,T|| )  [sec]  (2-113) 

(See  Equation  23  of  DNA-IR-82-022 . ) 

Next,  calculate  the  mean  weighted  distance  to  the  region  and  its 
standard  deviation: 

<Z>  =  sum(Z)/a2  [cm]  (2-114) 

o  <p 

sd(Z)  =  fsum(Z2)/a2  -  <Z>2]1/2  [cm]  (2-115) 

using  the  sums  calculated  in  Equations  2-97  and  2-98. 

The  mean  square  log  amplitude  fluctuation  (MSLAF)  is  required  in 
order  to  calculate  the  portion  of  the  total  phase  variance  whose  amplitude 
is  Rayleigh  distributed.  Calculate  the  MSLAF  using  the  primitive  inte¬ 
grals  of  Equations  2-43  through  2-45: 

x2  =  sum(Xo)  +  sum(Xl)  +  sum(X2)  (2-116) 

(See  Equation  A-3  of  DNA-IR-82-022 . )  Then  to  determine  the  portion  of  the 
phase  variance  which  is  Rayleigh  distributed,  let 

[c2(6_2n)J2-cS6"2n)  JT l/(6-2n)  n*3 

J,  =  (2-117) 

J2  +  In  ( c2 /c j )  n=3 
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where  J,  and  are  the  integrals  computed  in  Equations  2-53  and  2-54,  and 
cy  and  c,  are  constants:  ct  =  0.5,  c2  =  0.84.  (See  Eouations  A-9a  and 
A-9b  of  DNA-IR-82-02 2. )  And  let 

xs2  =  16(n-l)f(n)  c2n  J3  +  sum(Xl)  (2-118) 


where  sum(xi)  is  the  integral  computed  in  Eauation  2-44,  and  f(n)  is  the 
phase  weighted  mean  of  the  function  f(n,z)  defined  in  Equation  2-46: 


f(n)  - 


N 


d  a 


>  f(n,z.) 
f  1  dz 


A z  ]/ 


(2-119) 


The  phase  spectrum  is  arbitrarily  divided  into  two  regions  so  that 
fluctuations  whose  wave  number  is  smaller  than  some  critical  value  can  be 
considered  to  contribute  negligibly  to  the  total  phase  variance.  Denote 
the  value  of  the  MSLAF  integrated  up  to  that  wave  number  as  •>£•  The 
Rayleigh  threshold  is  usually  defined  to  be  x2  =  0.10.  PRPSIM  uses  the 
following  procedure  to  compute  the  Rayleigh  distributed  contribution  to 
the  phase  variance: 

If  x2  <  x2  or  (J,  <  0  and  a2/2  <  xf)  then 

a2  =  0  (2-120) 

else  if  sum(  xi)  £  x^  or  (J3  £  0  and  a^/2  _>  ^ )  then 

aR2  =  a2  [rad2]  (2-121) 

(where  sum(  xi)  is  the  integral  computed  in  Equation  2-44  and  o2  is  the 

v 

total  phase  variance  computed  in  Equation  2-40);  else  calculate 
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min  [ 


2J 


2 


T  •  °2J 


T  ?  =  < 


4n~l  a  (2n-2 )  <t>' 

c 


min  [ 


20i 


4n-l  ac(2n-2) 


-  ,  a2  ] 

\  *  a,  J 


n*3 


n=3 


where,  if  xc  <  x^ 


r 


{  (6-2n)  [. 


2  2 

Xc  -  xs 


a  =  < 
c  ' 


16c?n(n-l)f(n)J, 


xc  "  xs 

c2  exp  [ - - - ] 


-]  +  cf"2n} 


(6-2n ) 


n*3 


n=3 


16ct  (n-l)f (n)J1 


or,  if  X2  >  x2 
c  —  *s 


r 

{ (2-2n) 


V 


f (2-2n) 


2  2 
XC'XS 


2n 


16cinc2(n-l)f(n)J, 


-1  + 


(2-2n) 


n*3 


2  2 
Vxs 


16c?nc2(n-l)f(n)J1 


+  c 


2-2n, (2~2n) 


n=3 


(See  Equations  A-lla  through  A-13b  of  DNA-IR-82-022 . ) 


(2-122) 


(2-123) 


(2-124) 
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The  Su  scintillation  index  is  also  a  unitless  measure  of  the 
severity  of  the  amplitude  scintillation.  We  approximate  it  from  the  MSLAF 
by  usino  tnn  empirically  derived  relation 

Su  '  minimum  ( 1 .42( y2) 44,  1)  (2-125) 

where  ,  is  the  MTLAF  from  Equation  2-116.  (See  Equation  10  of  Reference 
h.  Nnt.e  that  this  approx  imat  ion  is  most  accurate  at  n  =  2. 

The  frequency  selective  bandwidth  is  computed  as 

U 

f,  =  f  r-2  +  —  ( r+e)  '-1/2  [Hz]  (2-126) 

c  R  K? 


where 


and 


f  is  the  signal  carrier  frequency  [Hz] 

K  is  the  signal  wave-number  [cm-1] 

-'n  is  the  Rayleiqh  phase  variance  from  Equations 
5  2-120  -  2-1??  [rad7] 

r  is  the  correction  factor  from  Equation  2-81  [cm-2] 

-  is  the  primitive  sum  from  Equation  2-78  [cm-2] 

Hm  =  ' sum(h)  1A/m-2  (2-127) 


is  the  finite  layer  correction  factor  derived  from  the  primitive  suin  from 
Equation  2-41.  (See  Equations  31-33  of  DNA- 1 R-82-02 2 . ) 
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The  Rayleigh  frequency  (fR)  is  defined  as  the  smallest  frequency 
of  the  plasma  oscillations  which  contributes  to  the  Rayleigh  phase  vari¬ 
ance: 


R  R 

f  P  ( f )df  =  2  r  p  (f)df  -  1  -  alio7  (2-128) 

■_f  *  o  5  R  * 

R 

where  ,72  is  the  Rayleigh  phase  variance,  a2  is  the  total  phase  variance, 

R  '  <J> 

and  Pt(f)  is  the  probability  distribution  in  frequency  for  the  electron 
density  fluctuations: 

P  (f)  =  2/rT  t  a2n'2  T  ( n.~  1 12  )  f  a2  +  (2m  f ) 2 1"  (n-1/2 )  (2-129) 

t  o  r(n-l)  0 

where  i0  is  the  signal  decorrelation  time  (Equation  2-113)  and  a  is  a 
unitless  coefficient  derived  from  the  phase  weighted  sum  of  the  Bn 
coefficient 

a  -  1  sum( h) I- (2-130) 

(See  Equations  2-29,  2-41,  and  Equation  39  of  DNA-IR-82-022  . )  For  con¬ 
venience,  we  define  the  [unitless]  scaled  Rayleigh  frequency 

ZR  =  2tt  (r0/a)fR  (2-131) 

Analytic  solutions  for  zR  exist  only  for  n  =  3/2  and  n  =  2.  For  n  =  2 

zR  =  !  rl  -  o2/a2l-2  -  1} -1/2  (2-132) 

For  other  values  of  n,  it  is  easy  to  show  that  the  limiting  solutions  for 
z  r  are 
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« 


1  im 

ZR  *  °>  ZR 


~  j/n  r(n-l) 

2  r (n-1/2) 


[1 


(2-133) 


and 


1  im 

ZR  +  00  *  ?R 


{ ( 2n-2) 


v'tt  r(n-l) 

2  r(n-l/2) 


(“R/o*>t'1/(2n'2> 


(2-134) 


Note  that  Equation  2-133  always  underestimates  the  true  value  of 
zR  whereas  Equation  2-134  always  overest imates  it.  To  compute  the 
true  value,  we  numerically  calculate  the  value  of  the  integral  (using 

Equation  2-128)  to  the  low  estimate  and  to  the  hiqh  estimate  of  zD, 

and  compare  these  quantities  to  the  known  value  of  the  integral.  We 
interpolate  logar ithmical ly  between  the  low  and  high  estimates  to  obtain  a 
new  estimate  of  zR,  and  numerically  calculate  the  integral  to  this 
value.  If  that  integral  exceeds  the  known  integral  to  the  true  value,  we 
replace  the  hiqh  estimate  of  zR  by  the  new  estimate  and  repeat  the 

process.  Otherwise  if  the  integral  is  smaller  than  the  known  integral,  we 
replace  the  low  estimate  of  z  by  the  new  estimate  and  repeat  the 

process.  In  practice,  2  to  4  iterations  are  usually  sufficient  to 

calculate  zD  to  within  1  percent  accuracy. 

K 


The  final  quantities  to  be  calculated  by  the  module  RFPROP  are 
the  means  and  standard  deviations  of  a  set  of  eiqht  parameters  which  are 
proportional  to  the  mean  and  standard  deviation  of  the  TEC  and  its  time 
derivatives.  The  variance  of  the  itln  order  time  derivative  of  the  TEC 
is  obtained  from  the  expression 


=  £  (t£_)2  a-™  /R 

"  cre  o 


P  (f)  (2nf)2j  df 

9 


(2-135) 
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where  P,(f)  and  a  are  defined  in  Equations  2-129  and  2-130,  and  f  is  the 

fp 

carrier  frequency  [Hz]  (see  Equation  A-21b  of  0NA-5662D8).  Equation  2-135 
can  be  expressed  more  conveniently  as 

c2  =  I  (!c_)2  a-m  (a_)2j  f  (z  )  (2-136) 

Nj  11  cre  To  ,J 

where 

f  .(z  )  =  /R  dz  z2j  (l+z2)-(n-l/2)  (2-137) 

n,J  K  0 

The  zeroth  order  derivative  may  be  evaluated  analytically  for  all  values 
of  n,  i .e. , 

a2  «  1  (!c_)2  a-m  ri-ag/o?]  (2-138) 

it  cr  K  $ 

e 

Analytic  solutions  for  the  higher  order  derivatives  exist  only  for  n  =  2 


f  2  ,  l  ( ZR )  =  109  ^  ZR  +  /1+zr)  '  d  +  zR2  )-  172  (2-139) 

fl,2  0R)  1/2  -  3f2>l(zR,!  (2-140) 

f2,,<2R>  =  |  !zJ(2-52R2)(lnR2)->/2  ♦  15f2tl(zR,]  (2-141) 


For  the  limiting  solutions  when  n  =  2 


1  im 

zR  ■*  0,  fp^. 


=  _1 _ 

<ZR)  (2 j+1 ) 


z2in 


j  =  1, 


2, 


(2-142) 
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1  im 

ZR  +  "*  f2,J'(zR) 

For  n  t  2,  no  analytic  solutions  exist  for  the  integral  Equation  2-137. 
PRPSIM  resorts  to  Simpson's  rule  to  solve  the  integral  numerically. 

The  mean  and  standard  deviataion  of  the  eight  additional  signal 
parameters  are  derived  from  the  mean  and  standard  deviation  of  the  TEC  and 
its  time  derivatives  as  listed  below: 

Mean  signal  phase  shift  and  its  standard  deviation: 

<phase  shift>  =  ACTEO  [rad]  (2-144) 

a(phase  shift)  =  A  a..  [rad]  (2-145) 

IN  n 

Mean  Doppler  shift  and  its  standard  deviation: 

<Doppl er  sh if t>  =  <TEC>  [Hz]  (2-146) 

2f r  dt 

^(Dopier  shift)  =  (~)  [Hz]  (2-147) 

Mean  Doppler  rate  and  its  standard  deviation: 


<Doppler 

rate> 

-  (—)  _ _  <TEC> 

2"  dt? 

[Hz  sec-1] 

(2-148) 

n{ Doppler 

rate) 

II 

Q 

TO 

[Hz  sec'1] 

(2-149) 

log(2zR)-l  j  =  1 


Mean  phase  "jerk"  and  its  standard  deviation: 


A  d 5 

'phase  jerk>  =  (  1  <TEC> 

2ir  dt3 

[Hz  sec-2] 

( 2-150) 

o(phase  jerk)  =  aN3 

[Hz  sec-2] 

(2-151) 

1  time  delay  and  its  standard 

deviat ion : 

rtime  delay>  =  B  <TEC> 

[  usee] 

(2-152) 

o(time  delay)  =  B  a^o 

[ usee] 

(2-153) 

Mean  time  delay  rate  and  its  standard  deviation: 

<t ime  delay  rate>  =  B—  <TEC>  [usee  sec'1]  (2-154) 

dt 

retime  delay  rate)  =  B  [usee  sec-1]  (2-155) 

Mean  time  delay  acceleration  and  its  standard  deviation: 

2 

/time  delay  acceleration>  =  B  — -  <TEC>  [gsec  sec-2]  (2-156) 

dt2 

a ( time  delay  accelerat ion)  =  B  a^2  [gsec  sec-2]  (2-157) 

Mean  time  delay  "jerk"  and  its  standard  deviation: 

•'time  delav  jerk>  =  B  —  <TEC>  [gsec  sec-  ]  (2-158) 

dt 3 

n(time  delay  jerk)  =  B  ^  [gsec  sec-3]  (2-159) 
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with 


A 


reX 


[cm2]  (2-160) 


B  = 


1x10 


cr 

6  e 


2*f; 


[cm2  gsec]  (2-161) 


where  c  is  the  speed  of  light,  rp  is  the  classical  radius  of  the 
electron,  f  is  the  carrier  frequency  [Hz],  and  X  is  the  wavelength  [cm]. 
This  point  marks  the  end  of  the  loop  in  frequency  in  which  the  signal 
parameters  are  calculated. 

Two  additional  subroutines  (RFTEMP  and  REFRAC)  are  called  by 
RFPROP  to  calculate  the  effective  temperature  of  the  receiving  antenna  due 
to  nuclear  enhanced  sky  temperature  and  to  calculate  the  bending  of  the 
ray  path  by  smooth  gradient  or  "gross"  refraction.  These  subroutines  will 
be  described  in  Sections  2.4  and  2.5  respectively. 

Table  3  summarizes  the  major  parameters  calculated  by  RFPROP  and 
indexes  the  numbers  of  the  equations  (in  Section  2-2)  pertinent  to  their 
calcul ation . 


44 


Table  3.  Summary  of  equation  numbers  for  major  parameters 
calculated  by  RFPROP. 


PARAMETERS 

UNITS 

EQUATION(s) 

decorrelation  distance 

[cm] 

106 

105 

decorrel ation  time 

[sec] 

113 

111 

109 

55 

Faraday  rotation  angle 

[rad] 

96 

frequency  selective  bandwidth 

[Hz] 

126 

78 

mean  square  log  amplitude 
fluctuation  (MSLAF) 

116 

45 

44 

43 

mean  signal  phase  shift 

[rad] 

144 

mean  signal  time  delay 

[usee] 

152 

mean  square  angle  of 
energy  arrival 

[rad2] 

108 

107 

jjiean  weighted  distance 

[cm] 

114 

97 

Rayleigh  phase  variance 

[rad2] 

122 

121 

120 

Rayleigh  frequency 

[Hz] 

128 

S4  scintillation  index 

125 

total  absorption 

[dB] 

93 

90 

total  electron  content  (TEC) 
and  time  derivatives 

[cm-2  seed] 

94 

23 

22 

19 

total  electron  standard  devia¬ 
tions  and  tire  derivatives 

[cm-2  sec_J] 

136 

total  phase  variance 

[rad2] 

40 

38 

45 


2.3 


ANTENNA  FILTERING. 


The  subroutine  ANTMOD  is  called  by  PRPSIM  to  calculate  the  modi¬ 
fication  of  the  signal  parameters  resultinq  from  the  loss  of  signal  power 
caused  by  angular  scattering  at  the  receiving  and  transmitting  antennas. 
The  equations  used  in  ANTMOD  are  taken  from  those  derived  by  Dana  (1985). 9 
The  derivations  assume  that  both  antennas  are  character ized  by  Gaussian 
beam  profiles  and  that  both  antennas  are  pointed  along  the  line  of  sight. 
Four  arrays  of  beamwidths  are  passed  to  ANTMOD  by  RFPROP .  At  each  propaga¬ 
tion  freauency,  each  antenna  is  character  ized  by  two  beamwidths  along  two 
orthogonal  axes  perpend icu 1 ar  to  the  axis  of  the  antenna.  The  beamwidth 
is  defined  as  the  anqle  [rad]  subtended  by  the  half  power  (3  dB)  points  on 
the  directive  gain  function  on  opposite  sides  of  the  antenna  axis.  Cur¬ 
rently,  only  circularly  symmetric  antennas  are  allowed  by  RFPROP,  and  the 
beamwidths  in  the  two  orthoqonal  directions  are  equal.  Here  we  denote  the 
beamwidths  nf  the  transmitting  and  receiving  antennas  along  the  X  and 
Y-ayes,  respectively,  as  °Tx ,  ,  and  .  The  signal  loss  caused 

bv  scattering  is  given  by 

1  =  10  looiobx)  +  10  log !  o  (Ay)  [dB]  (2-162) 

whore 


Ax)‘  =  1  +  k' 


■y  / 


l  +  k  1 


/’Tx  + 

nRx 

rv 

Tx 

Rx 

0 

? 

r~Ty  + 

'«y 

-r  + 

°Ty 

Py 

(2-163) 


(2-164) 


when.  -1  ,,nd  V  are  the  mean  square  angles  of  energy  arrival  in  the  x 
Tx  Tv  ,  , 

,.j n . !  y  <\  irer  t.  ir.ns  at  the  transn. i  t ter  ,  and  ip  are  the  analogous  quanti¬ 
fy  j«.-,  af  til,,  rc-ro  i/er,  and  k  '  is  the  constant  Pln(2)  =  5.645.  For  a 
c  ir^ujar  ant'  nna  w  r. an  replace  the  mean  square  angles  in  the  x  and  y 
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directions  by  the  minimum  and  maximum  values  in  the  p  and  q  directions 

calculated  in  Equations  2-107  and  2-108.  A  beamwidth  of  it  (180°)  or 

qreater  is  interpreted  by  ANTMOD  as  an  omn i -d irect iona 1  antenna  with  zero 

scattering  loss  ((Ax)2  =  (Ay)2  =  1).  ANTMOD  also  computes  the  total 

signal  loss  which  is  simply  the  sum  of  the  absorption  and  the  scattering 
loss.  The  apparant  frequency  selective  bandwidth  and  decorre 1  at  ion  time 
are  increased  by  scattering: 


=  Ax  Ay  f0 

[Hz] 

(2-165) 

=  Ax  Ay  tq 

[sec] 

(2-166) 

where  f 0  and  r0  are  the  antenna  independent  values  calculated  by  Equations 
2-126  and  2-113,  respectively.  Note  that  the  degree  of  antenna  filtering 
is  frequency  dependent. 

2.4  ANTENNA  NOISE  TEMPERATURE. 

PPPSIM  calculates  the  effective  antenna  noise  temperature  of  the 
receiving  antenna  due  to  a  hot,  ionized  medium.  The  calculation  includes 
only  the  contribution  of  the  medium.  All  other  sources  (solar,  cosmic, 
the  earth,  etc.)  are  not  considered  here.  In  theory  the  temperature  must 
be  integrated  over  every  point  in  the  grid  weighted  by  the  antenna  gain 
function  in  the  direction  of  the  point.  Such  an  analysis  is  beyond  the 
scope  of  a  propagation  code  whose  primary  purpose  is  the  calculation  of 
scintillation  effects.  We  compromise  by  calculating  the  temperature  with 
the  "beam-f i 1 1 inq"  approx imat ion  formula  at  five  points  on  the  gain  func¬ 
tion,  and  computinq  the  mean  and  standard  deviation  of  these  five  values. 
One  direction  is  along  the  line  of  sight  (LOS)  directed  from  the  receiver 
to  the  transmitter.  The  other  four  directions  are  displaced  one-half 
beamwidth  from  the  LOS  along  two  orthogonal  directions  perpend icu 1 ar  to 
the  LOS.  As  before,  the  beamwidth  of  the  antenna  is  defined  to  be  the 
ariqle  subtended  by  the  half  power  (3  dB)  points  on  the  directive  gain 
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function  on  opposite  sides  of  the  antenna  axis.  Currently,  we  assume  that 
the  antenna  gain  function  is  circularly  symmetric  and  that  the  antenna 
axis  is  oriented  along  the  LOS.  RFPROP  may  easily  be  generalized  to 
accomodate  non-circular  antennas,  however.  A  warning  flag  is  printed  if 
the  mean  does  not  exceed  the  standard  deviation  by  a  factor  of  two  or 
greater.  The  flag  is  intended  to  indicate  that  the  variation  of  the 
temperature  across  the  antenna  gain  function  may  be  so  large  that  the 
small  number  of  points  used  may  yield  an  inaccurate  result.  Although  the 
absence  of  a  flag  indicates  a  reliable  temperature  calculation,  a  flag 
does  not  necessarily  indicate  that  the  resu’t  is  inaccurate.  Even  when 
the  result  is  flagged,  the  mean  is  a  much  better  estimate  of  the  true 
antenna  temperature  than  the  single  beam-filled  approximation  calculated 
along  the  LOS  only. 

The  antenna  temperature  calculation  is  performed  by  a  submodule 
of  RFPROP  called  RFTEMP .  RFTEMP  is  called  once  by  RFPROP  for  each  link 
simulation.  The  arguments  passed  are  the  vectors  containing  the  coordi¬ 
nates  of  the  transmitter  (path  origin)  and  the  receiver  (path  end),  the 
number  of  RF  carrier  frequencies ,  the  array  containing  those  frequencies, 
the  array  containing  the  effective  beamwidths  along  the  X-axis  at  those 
frequencies,  the  array  of  beamwidths  along  the  Y-axis  (currently  these 
equal  those  along  the  X-axis),  and  two  integers  which  specify  the 
sequential  indices  of  the  first  and  last  point  (relative  to  the 
transmitter)  whose  mean  local  electron  density  or  standard  deviation 
exceeds  an  arbitrary  threshold.  (Currently  the  threshold  is  1  cm'  for 
both  parameters.) 

The  antenna  temperature  calculation  can  be  the  most  time 
consuming  of  all  PRPSIM  calculations,  and  for  that  reason  it  is  optional. 
Data  retrieval,  rather  than  the  complexity  of  the  calculations, 
constitutes  the  major  demand  for  CPU  time  .  Each  of  the  five  paths 
requires  approximately  the  same  number  of  data  retrievals  required  for  all 


the  computuat ions  performed  by  RFPROP.  Furthermore,  because  the  antenna 
qain  function  and  beamwidth  are  functions  of  frequency,  a  set  of  four 
different  auxiliary  paths  must  be  defined  for  each  operating  freauency. 
If  n  operating  frequencies  are  specified  for  a  link,  it  is  easy  to  show 
that  the  number  of  data  retrievals  will  increase  by  a  factor  of  4n+?  if 
antenna  temperature  is  calculated.  (In  theory,  this  ratio  could  be 
reduced  to  4n+l  if  the  temperature  along  the  central  path  were  to  be 
calculated  by  RFPROP.  Owing  to  the  directional  asymmetry  of  the  calcula¬ 
tion  of  temoerature  and  several  of  the  other  parameters,  it  is  simpler  for 
numerical  reasons  to  calculate  the  temperature  along  the  central  path  in 
RFTEMP.)  We  do  reduce  the  CPU  time  somewhat  by  eliminating  the  construc¬ 
tion  of  the  auxiliary  paths  for  situations  in  which  it  is  easily  demon¬ 
strated  that  the  antenna  beamwidth  is  filled  by  a  homogeneous  ionized 
region. 

RpJFMP  calls  PATHPT  to  qenerate  a  set  of  points  along  the  propagation  path 
from  "infinity"  (at  a  point  on  the  path  located  well  outside  the  phenome¬ 
nology  grid)  to  the  receiver.  The  indices  (relative  to  the  receiver)  of 

o 

the  first  and  last  point  whose  mean  density  exceeds  the  1  cm'  threshold 
are  also  computed.  The  code  then  begins  to  loop  over  those  points  alonq 
the  path  whose  indices  lie  between  the  threshold  indices  to  calculate  the 
incremental  contribution  to  the  total  antenna  temperature  from  the  plasma 
in  that  reqion.  The  submodule  DATAPT  (described  in  Section  3.4)  is  called 
to  retrieve  the  set  of  interpolated  environmental  data  for  each  point. 
The  data  set  returned  is  truncated  to  only  those  data  needed  for  the 
temperature  calculation:  (1)  the  mean  electron  density  [cm-3],  (2)  the 
electron  density  standard  deviation  Tcm'3],  (3)  the  rms  electron  density 
[cm""],  (4)  the  electron  temoerature  [°K],  (5)  the  electron-neutral  colli¬ 
sion  frequency  [sec-1],  and  (6)  the  array  of  electron-ion  collision 
frequencies  (sec'1].  (See  Table  4  of  Section  3.4.) 
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If  we  designate  the  local  ion  temperature  of  a  plasma  layer  of 
thickness  Az-j  as  ti,  it  can  be  shown  that  the  temperature  of  the 
receiving  antenna  is 

dTi  =  t.  n  -  exp(-2Kv.  Az  • )  lLpi  (2-167) 

or 

dT  .  =  t.  4  -  10"(aiAZl'/10)lLpi  (2-168) 

if  the  layer  is  of  uniform  temperature  and  completely  subtends  the  beam- 
width  of  the  antenna.  (See  Equation  27  of  MRC-R-156,  Volume  II7.)  Lp  •  is 
the  loss  caused  by  the  absorption  of  the  intervening  material  between  the 
region  and  the  receiver: 

i 

Lp.  =  exp  1-2  ^  Kv.Az.l  (2-169) 

j  J  J 

or 

i 

-r  l  ntjAZj/10’1 

Lpi  =  10  J  '  (2-170) 

where  Kv.  is  the  linear  coefficient  of  absorption  [cm-1],  and  a  is  the 
coefficient  of  absorption  which  was  derived  in  Equation  2-93  of  Section 
2.2.  Kv  and  n  are  related  by 

't  =  — -  Kv  [dB  cm-!]  (2-171) 

In  (10) 

The  factor  of  two  enters  the  exponential  arguments  because 
temperature  is  a  measure  of  power  whereas  the  coefficient  of  absorption  is 
defined  in  terms  of  amplitude  attenuation  as  a  function  of  distance.  Note 
that  Kv  is  a  function  of  both  frequency  and  location.  Thus,  it  is  neces¬ 
sary  to  calculate  the  temperature  increment  for  each  carrier  frequency  and 
at  each  point  along  the  path.  The  temperature  of  the  receiver  antenna  (as 
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a  function  of  frequency)  is  calculated  by  summing  over  all  points  along 
the  path 
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(2-172) 


Note  that  Equation  2-172  depends  on  the  order  of  summation.  Because  we 

want  the  temperature  at  the  receiver,  we  must  sum  from  the  receiver  (i=l) 

toward  "infinity"  ( i=N) .  We  also  calculate  the  mean  distance  to  the 
reqion  weiqhted  by  the  incremental  temperature: 

N  N 

<  z>  j  =  y  z-  dT  ■/  dT  .  [cm]  (2-173) 

Equation  2-172  is  correct  only  if  the  medium  is  uniform  over  the  beamwidth 
of  the  antenna.  To  test  whether  this  is  the  case,  we  assume  that  the 

SCENARIO  grid  resolution  is  sufficiently  good  to  assure  that  the  proper¬ 
ties  of  the  medium  are  reasonably  constant  across  a  grid  cell;  and  if  we 
can  show  that  the  grid  spacing  at  the  most  highly  ionized  region  along  the 
LOS  is  greater  than  the  arc  subtended  by  the  beamwidth  at  that  distance, 
then  the  beam-filling  approximation  may  be  used  in  lieu  of  integration 
over  the  gain  function.  First,  using  a  simple  search  algorithm  and  the 

array  of  grid  coordinates  read  and  stored  from  the  SCENARIO  data  file,  we 

calculate  the  indices  of  the  cell  containing  the  point  at  the  temperature 
weighted  distance  to  the  region  computed  in  Equation  2-173.  The 
dimensions  of  the  cell  along  the  three  orthogonal  coordinate  axes  are 

computed  as 

dx  .  =  h i  du  •  [cm]  (2-174) 
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where  du •  are  the  differences  between  curvilinear  coordinates  read 
from  the  SCENARIO  file,  and  h-  are  the  scale  factors  as  calculated 
from  Equations  3-40,  3-41,  and  3-42  of  Section  3.4.  The  angle  subtended 
by  the  smallest  of  these  dimensions  at  the  temperature  weighted  distance 
to  the  cell  is 


e 


rmin(dx1  ,dx2,dx3) 

2  arctan  i - - - - - _l 

2<z>t 


(2-175) 


If  this  anqle  is  greater  than  the  beamwidth  of  the  antenna,  the 
ionized  medium  is  considered  homogeneous  throughout  the  beamwidth  of  the 
antenna  so  that  the  beam-filling  approximation  is  valid.  In  that  case  we 
do  not  construct  the  four  auxiliary  paths  as  outlined  below,  and  the 
antenna  temperature  is  equated  to  the  value  obtained  with  Equation  2-172. 
The  standard  deviation  of  the  temperature  across  the  antenna  face  is  set 
at  zero. 


RFTEMP  then  loops  through  the  set  of  carrier  frequencies.  If  e 
does  not  exceed  the  antenna  beamwidth  for  a  given  frequency,  we  proceed  to 
construct  a  set  of  four  auxiliary  paths  for  that  frequency.  The  paths  are 
directed  toward  the  half-power  contour  of  the  gain  function  and  spaced  at 
ninety  degree  increments  in  azimuth  about  the  LOS  (which  is  assumed  to  be 
coincident  with  the  antenna  normal).  Define  the  path  elevations  and 
azimuths  relative  to  the  LOS: 

tt/2  -  1  beamwidth  i  =  1  or  3 

t.  =J  2  *  (2-176) 

tt/2  -  —  beamwidth  i  =  2  or  4 

2  y 


and 
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(2-177) 


<t> ■  =  0,  tt/2,  it,  3 tt/ 2  for  i  =  1,2, 3, 4 

(Note  that  the  elevations  of  the  auxiliary  paths  relative  to  the  LOS  are 
frequency  dependent.)  The  end  of  each  auxiliary  path  is  located  at 
"infinity"  (i.e.,  outside  the  grid).  All  computations  are  done  in  the 
tangent  plane  (east,  north,  vertical)  system  of  the  antenna.  To  transform 
the  coordinates  of  the  four  auxiliary  paths  from  antenna  coordinates  to 
tangent  plane  coordinates:  (1)  rotate  the  system  counterclockwise  about 
the  W-axis  (normal  to  the  antenna  face)  through  the  angle  $  to  align  the 

V-axis  in  the  horizontal  plane,  (2)  rotate  counterclockwise  about  the 
V-axis  through  the  angle  n/2-e  to  align  the  W-axis  with  the  Z-axis 
(vertical)  and  (3)  rotate  counterclockwise  about  the  Z-axis  through  the 
angle  tt/2+X  to  align  the  U-axis  with  the  X-axis  (east)  and  the  V-axis  with 
the  Y-axis  (north).  (See  Figure  3.)  Currently,  only  circular  antennas 
are  allowed  by  RFPROP,  and  i>  is  undefined;  for  non-circular  antennas,  the 
U-axis  woul'4  -'^rmally  be  defined  to  lie  along  the  major  axis  of  the  gain 

function,  and  >| >  would  be  defined  as  the  angle  (measured  clockwise)  from 

the  horizontal  to  the  V-axis. 

RFTEMP  tests  to  determine  whether  any  of  the  auxiliary  paths 
intersects  the  earth.  If  a  path  does  intersect  the  earth,  its  length  is 

redefined  to  be  the  distance  from  the  receiver  to  the  (nearest) 
intersection  point.  For  each  of  the  four  auxiliary  paths,  let 

b  =  i  rz0  +  (r2-r2)/z0l  (2-178) 

2  * 

and 

c  =  r2  -  (2-179) 

where  z0  is  the  original  length  of  the  auxiliary  path,  r  is  the  radial 
distance  to  the  receiver  from  the  center  of  the  earth,  r^  is  the  radial 


53 


Figure  3 


Geometric  transformation  from  line  of  sight  (LOS)  coordinates 
to  antenna  tangent  plane  (east,  north,  vertical)  coordinates. 
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distance  to  the  transm i tter ,  and  rg  is  the  radius  of  the  earth.  If  b  <  0 
or  if  c  >  b2,  the  path  does  not  intersect  the  earth,  and  we  proceed  to 
calculate  the  temperature  by  integrating  along  the  entire  length  of  the 
path.  Otherwise,  calculate  the  distance 


x  =  b  -  /b2-c  (2-180) 

If  x  <  z0,  the  path  intersects  the  earth,  and  we  redefine  the  integration 
path  lenqth  to  be  x. 

PATHPT  (see  Section  3.3)  is  then  called  to  generate  a  set  of 
grid  intersection  points  for  the  auxiliary  path.  The  incremental 
temperature  is  calculated  and  summed  along  the  auxiliary  path  in  exactly 
the  same  manner  as  along  the  LOS  except  that  there  in  no  inner  loop  in 
frequency  for  the  path  point.  Frequency  is  the  outermost  loop  for  the 
auxiliary  paths  rather  than  the  innermost.  The  mean  antenna  temperature 
and  the  standard  deviation  are  calculated  by  using  the  (frequency 

dependent)  temperature  integrated  along  the  LOS  and  the  four  auxiliary 

paths : 

4  4 

<T:>  -  V  W ,T  -/  V  w.  [°K]  (2-181) 

i=0  i=0 

4  4 

aT  =  f  V  w.T2/  V  W.  -  <T>2]1/2  [°K]  (2-182) 

1  i=0  1  1=0 

where  T0  is  the  temperature  calculated  along  the  LOS  for  that  frequency, 
T[  through  T,+  are  the  temperatures  calculated  along  the  auxiliary  paths 
and  W-  are  the  weiqhting  parameters.  Thp  gain  function  along  the 
auxiliary  paths  is  half  that  along  the  LOS,  implying  that  we  might  choose 

W0  =  1  and  W-  =  0.5  for  i  =  l  to  4.  The  contribution  of  each  of  the 

auxiliary  paths  to  the  total  solid  angle  is  somewhat  greater  than  that  of 
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the  LOS,  however,  so  that  the  kL  should  be  somewhat  larger  than  0.5  (but 
smaller  than  1).  Here,  we  have  conservatively  set  all  the  weighting 
parameters  to  1,  even  though  this  choice  probably  biases  the  estimate  of 
the  standard  deviation  toward  values  somwhat  greater  than  its  real  value. 

2.5  REFRACTION. 

PR  PS  I  rd  calculates  the  refractive  effects  on  an  RF  signal  propa¬ 
gating  in  an  i nhomogeneous  ionized  medium  whose  properties  are  known  only 
at  specified  points.  The  primary  objective  is  to  provide  reasonably 
accurate  calculations  of  refraction,  multipath,  focusing  and  defocusing, 
phase  shifts,  and  group  delay  in  the  disturbed  environment  modelled  by 
SCENARIO.  The  refractive  effects  considered  here  are  those  effects  caused 
by  large-scale  ionization  structure  in  the  propagation  medium  -  often 
referred  to  as  "gross"  or  "smooth  gradient"  refractive  effects.  That  is, 
the  index  of  refraction  is  a  sufficiently  slowly-varying  function  in  space 
so  that  geometrical  optics  or  "ray  theory"  is  applicable;  and  furthermore, 
the  index  of  refraction  is  a  reasonably  smooth  and  well-behaved  function 
in  the  regions  through  which  most  of  the  received  signal  energy  propa¬ 
gates.  The  approximation  of  geometrical  optics  imposes  no  significant 
constraint  because,  when  the  refractive  index  is  well-behaved,  this 
approximation  requires  that  the  properties  of  the  medium  change  very 
little  over  dimensions  comparable  to  a  wavelength.  At  frequencies  in  the 
UHF  band  and  higher,  this  requirement  is  easily  satisfied  whenever  signals 
can  propagate  at  all.  Section  2.2  treats  the  effects  of  random  inhomo- 
geneities  in  the  index  of  refraction  due  to  relatively  small-scale  ioniza¬ 
tion  structure  encountered  in  the  striated  plasma. 

2.5.1  Primary  Angle  Bending  and  Multipath. 

The  technique  used  here  (implemented  in  routine  REFRAC  and  its 
associated  subroutines)  to  compute  refractive  effects  on  RF  signal  propa¬ 
gation  is  one  derived  in  Reference  10.  The  technique,  referred  to  as  the 
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multiscreen  method,  is  more  soundly  based  in  propagation  theory  than  the 
simpler  refraction  computational  methods  used  in  system,  analysis  codes, 
and  it  is  much  more  efficient  than  ray-tracinq  techniques. 

The  multiscreen  method  is  a  stationary  phase  approach  in  which 
the  phase  shifts  caused  by  the  propagation  medium  are  represented  by  a 
series  of  phase  screens  along  the  propagation  path.  Each  phase  screen  is 
constructed  by  means  of  a  Taylor  series  expansion  about  the  line  of  siqht 
(LOS)  at  each  point  at  which  environmental  data  are  available.  The 
geometry  of  the  primary  propagation  path  is  estimated  by  applying  Fermat's 
principle  to  a  simplified  eouation  for  the  phase  path  length  throuqh  the 
multiple  phase  screens.  The  simplification  is  such  that  the  condition  for 
stationary  phase  requires  the  solution  of  2n  simultaneous  linear  equa¬ 
tions,  where  n  is  the  number  of  phase  screens.  The  set  of  2n  equations  is 
a  five-diagonal  system  that  is  easily  solved  by  a  modification  of  the 
Grout  reduction  technique.  The  resulting  path  geometry  is  considered  a 
1st  order  approximation  to  the  actual  path  qeometry.  A  parameter  that 
describes  the  perturbation  of  the  actual  path  qeometry  from  the  represent¬ 
ative  geometry  is  then  introduced.  A  more  complete  equation  for  the  phase 
path  length  of  the  perturbed  path  through  the  multiple  phase  screen  is 
formulated  using  hiqher-order  terms  in  the  series  expansions.  The  final 
solution  for  the  stationary  paths  involves  a  cubic  equation  in  terms  of 
the  perturbation  parameter,  which  yields  up  to  three  possible  propagation 
paths. 

The  phase  path  length  of  a  ray  path  is  given  by  the  following 

i nteqral 

R 

3  =  -  -  ;  udr  [radl  (2-183) 

c  o 

where  j  is  the  anqular  signal  frequency  in  rad  sec-1,  c  is  the  speed  of 
light  (2.999F+10  cm  sec-1),  and  u  has  the  following  form: 
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2 

u  -  (1  -  ^H)1/2  (2-184) 

u2 

where  u  is  the  angular  plasma  frequency  in  rad  sec'1.  Fermat's  principle 
states  that  the  paths  along  which  the  signal  propagates  are  paths  of 
stationary  phase  -  that  is,  paths  for  which  the  phase  is  constant  for 
small  perturbations  of  the  path  geometry.  The  path  geometry  is  the 
unknown  to  be  determined.  Hence,  one  first  assumes  some  kind  of  reason¬ 
able  geometry  as  a  function  of  one  or  more  independent  parameters,  and 
then  solves  for  the  parameter  values  that  satisfy  the  condition  of  sta¬ 
tionary  phase. 

In  setting  up  the  path  geometry,  we  note  that  the  information 
about  the  propagation  medium  consists  of  the  mean  electron  density  and  its 
spatial  derivatives  (described  in  detail  in  section  3.4)  at  a  sequence  of 
points  along  the  LOS  path  (described  in  detail  in  section  3.3)  to  the 
receiver.  In  effect,  we  have  information  on  n  phase  screens  along  the 
path,  as  illustrated  in  figure  4.  The  points  r  ,  r  ,  ...,  rn  are  calcula¬ 
tion  points  along  the  LOS  to  the  receiver  located  at  rR£VR  •  The  electron 
density  and  its  transverse  spatial  derivatives  up  to  some  reasonable  order 
are  known  at  each  point.  A  possible  refracted  path  is  then  postulated. 
The  refracted  path  intersects  the  planes  normal  to  the  LOS  at  x.,  y ■  ,  i  = 
1,  2,  ...,  n.  The  x-axis  lies  in  the  "vertical"  plane  and  the  y-axis  lies 
in  the  plane  normal  to  the  "vertical"  plane  (traverse  plane).  (Note  that 
the  x  and  y  axes  here  are  the  same  as  the  U  and  V  axes  defined  in  section 
2.2.1.) 

With  this  geometry  the  phase  can  be  evaluated  from  equations 
(2-183)  and  (2-184).  Under  the  condition  that  w2  »  w2,  which  is  well 
satisfied  whenever  absorption  is  not  catastrophically  1 arge  at  the 
frequencies  of  interest,  and  assuming  small  refractive  bendinq  angles,  the 
total  phase  shift  along  the  refracted  path  is  given  by: 
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Figure  4.  Multiple  phase  screen  refraction  geometry. 


t  he  prop t.  i ori  p rj t  h s  must,  connor t.  t.hp  t.ransm  i  1 1  or  finrl  ror.oivor, 
wo  h.ivo  t.ho  following  bouruMry  corol  i  t  i  oris : 
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x0  =  yo  =  0 
rn+l  =  rRCVR 
xn+l  yn+l 


(2-188) 


The  electron  density  at  the  intersection  of  the  refracted  path 
with  each  screen  can  be  evaluated  by  means  of  a  Taylor  series  expansion. 
To  linearize  the  set  of  2n  simultaneous  equations,  only  terms  through 
second  order  are  retained  in  the  expansion  at  this  stage.  (A  method  of 
incorporating  the  effects  of  higher-order  terms  is  introduced  shortly.) 
Thus  the  system  of  2n  equations  to  be  solved  is: 
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where  the  subscript  i  denotes  values  calculated  at  the  point  r.  on  the 
LOS. 


If  these  2n  equations  and  the  unknowns  are  arranqed  in  the 

orders  x  ,  v  ,  x  ,  y  ,  ...,  x  .  y„  ,  the  coefficient  matrix  has  a  five- 
l’  l’  2  2  n’  Jn  ’ 

diagonal  form.  This  five-diagonal  system  can  be  readily  solved  using  the 
Crout  reduction  method  (see,  for  example,  reference  11).  The  system  of 
equations  can  be  written  as: 
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From  Equations  2-189,  the  elements  of  the  auqmented  matrix  of  the  system 
are  qiven  by: 
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It  can  be  readily  verified  that  the  Crout  reduction  preserves 
the  five-diagonal  property,  and  that  the  elements  of  the  auxiliary  matrix 
are  qiven  by  the  simplified  forms: 
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Within  the  limits  of  the  small-anqle  approximations  and  the 
neqlect  of  hiqher-order  derivatives  in  the  Taylor  series  expansions,  these 
solutions  for  x.  and  V-  provide  an  essentially  exact  solution  for  the 
propagation  path  when  an  adeouate  number  of  calculation  points  are  used. 


The  simplest  method  of  incorporat inq  the  effects  of  hiqher-order 
spatial  derivatives  in  the  Taylor  series  expansions  is  to  assume  that  the 
solution  just  obtained  is  typical  of  the  propagation  path  geometry,  and 
that  the  actual  path  (or  paths)  can  be  represented  by  a  perturbation  of 
this  geometry.  In  this  case  we  consider  the  actual  propagation  path  to 
intersect  each  phase  screen  at  the  point  (5x.,  6y.),  where  x^  and  y-  are 
the  values  just  obtained  from  Equations  2-191,  and  6  is  a  parameter  whose 
value  is  to  be  determined.  Since  only  a  sinqle  unknown  parameter  is  now 
involved,  terms  in  the  Taylor  series  expansions  through  fourth  order  can 
be  retained,  which  yields  a  sinqle  cubic  equation  in  6  to  solve  for  the 
stationary  phase  paths: 
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At  this  point  we  note  that  if  all  the  third  and  fourth  spatial 
derivatives  of  the  electron  concentration  were  zero,  then  s2  and  s3  would 
vanish  and  <5  would  be  equal  to  -Sq/s^  However,  under  this  condition  we 
know  that  5  must  equal  unity,  since  the  unknown  x.  and  y.  give  the  sta¬ 
tionary  phase  path  when  the  third  and  fourth  derivatives  are  neglected. 
Therefore,  we  have  the  additional  information  that  so  =  -si,  which  can  be 
used  to  simplify  the  computation  or  to  check  the  multiscreen  solution. 
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The  cubic  equation  (Equation  2-192)  can  be  solved  by  standard 
techniques.  The  values  of  the  perturbation  parameter  6  thus  obtained  are 
used  to  modify  the  original  multiscreen  solutions  to  obtain  improved 
solutions  for  the  propagation  paths.  The  procedure  can  yield  any  number 
of  propagation  paths  from  zero  to  three,  depending  upon  the  electron 
density  gradients.  Thus,  refractive  multipath  conditions  are  directly 
treated  by  this  procedure. 

If  one  is  to  consider  the  transmitter  as  a  radar  and  the  re¬ 
ceiver  as  a  target,  when  more  than  one  propagation  path  exists  between  the 
radar  and  tarqet,  say  m  paths,  then  m2  returns  are  present  at  the  radar 
receiver  resulting  from  each  combination  of  outqoing  and  incoming  paths, 
interference  between  the  multiple  returns  gives  rise  to  amplitude  fading 
and  jitter  in  the  radar  measurements. 

Usinq  the  values  of  <5  in  conjunction  with  the  unperturbed  screen 
solutions,  the  refractive  bendinq  angles  at  each  end  of  the  oropaqation 
path  are  given  by  (see  figure  4): 
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The  orientation  of  these  refractive  anqles,  relative  to  the  vertical 
direction  (x-axis),  measured  clockwise  as  seen  from  the  transmitter  end  of 
the  oath,  are  given  by: 
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2.5.2  Focusing  and  Defocusing. 

Focusing  and  defocusinq  can  also  occur  in  a  refractive  environ¬ 
ment.  In  such  cases  the  ionosphere  acts  somewhat  as  a  lens,  causing 
either  enhanced  or  reduced  signal  strength  at  the  receiver  depending  on 
the  character i st ics  of  the  spatial  variation  of  the  refractive  index. 
Focusing  and  defocusing  are  usually  significant  only  when  a  propagation 
mode  is  just  on  the  verge  of  disappearing  or  reappearing. 

Assume  that  the  first  path  solution  to  the  receiver  (located  at 
rn+l  =  rRCVR’  xn+l  =  yn+l  =  has  keen  °ktained,  and  the  coefficients 
of  the  cubic  perturbation  equation  have  been  computed.  To  estimate  the 
focusing/defocusing  properties  of  the  medium  in  the  refractive  gradient 
direction,  another  propagation  path  solution  is  obtained  for  a  point 
displaced  slightly  in  the  receiver  plane  in  the  direction  in  which  the 
path  intersection  with  the  last  screen  is  displaced. 


Thus,  for  the  gradient,  direction,  we  take  the  displaced  point 
to  be  located  in  the  receiver  plane  at 
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where  e  is  a  constant  displacement  from  the  receiver,  and  x  ,  yn  are  the 
coordinates  of  the  point  at  which  the  previously  computed  path  to  the 
receiver  intersected  the  last  screen.  Note  that  these  values  of  xn+^, 
yn+j  are  constants,  not  variables,  in  the  equations  to  be  solved  for  the 
propagation  paths  to  this  displaced  point. 


Similarly,  in  the  direction  normal  to  the  gradient  direction,  we 
take  the  displaced  point  in  the  receiver  plane  at 
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In  either  case,  the  resulting  system  of  simultaneous  linear 
equations  that  must  be  solved  to  obtain  the  propagation  path  to  the 
disp.aced  point  is  identical  to  Equation  2-189,  with  the  exception  that 
xn+l’  yn+l  dre  no  'on9er  zero-  fact,  the  only  changes  in  the  elements 
of  the  augmented  matrix  of  the  system  are  that  the  final  two  members  of 
the  column  vector  c  each  contain  one  additional  term: 


rn+l  -  Vl  3nen 


n+1 


'2n-l 


4N 


3x 


rn+l  ‘  rn 


(2-205) 


rn+l  ~  rn-l  3nen 

c~  = - -  + 

2n  4N  By 


'n+1 


r  -  r 
n+1  n 


(2-206) 


With  the  Grout  reduction  method  of  solution,  it  is  seen  that  the  above 
affects  only  the  last  two  values  of  c'  in  the  auxiliary  matrix: 
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These  two  changes  propagate  back  and  modify  all  of  the  x.. ,  y. ,  i=l,2, 
...,n.  The  focusing/defocusing  properties  in  each  of  the  two  directions 
are  then  computed  by  determining  the  angular  separations  between  the 
receiver  path  A’  J  the  displaced  paths  compared  to  those  which  would  result 
in  free-space. 

Once  the  screen  equations  are  solved  for  the  displaced  path,  the 
one-dimensional  one-way  defocusinq  ratio  in  the  refractive  qradient 
direction  is  computed  as: 
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where  primes  denote  values  obtained  for  the  displaced  focusing  path.  For 
the  direction  orthogonal  to  the  qradient  direction  a  similar  solution  is 
obtained.  (Note  that  a  defocusing  ratio  greater  than  one  implies  a  loss 
of  signal  enerqy  and  the  complement  of  this  condition  implies  a  focusing 
or  increase  in  signal  power  measured  at  the  receiver.) 

2.5.3  Refractive  Effects  on  Phase  and  Time  Delay. 

Refractive  bending  of  the  propagation  path  causes  the  phase  and 
group  path  lengths  to  be  changed  from  the  values  they  would  have  if  the 
wave  propagated  along  the  LOS  to  the  receiver.  These  chanqes  can  be 
readily  calculated  using  the  multiscreen  method. 
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When  the  multiscreen  solution  for  the  primary  propagation  path 
has  been  obtained  and  the  coefficients  of  the  cubic  perturbation  equation 
have  been  calculated,  the  phase  shift  along  each  of  the  possible 
propagation  paths  is  determined  as  follows.  Integration  of  Equation  2-192 
with  respect  to  6  yields  the  phase  path  length  for  each  of  the  paths. 
This  in  conjunction  with  Equation  2-183  gives  the  phase  shift  along  the 
refracted  paths: 

9  =  —  f s o 6  +  I  s ! 62  +  I  s263  +  I  s354  )  +  9q  [rad]  (2-210) 

where  90  is  the  phase  shift  evaluated  along  the  LOS  to  the  receiver, 
including  the  effects  of  electron  content  along  the  LOS.  When  absorption 
is  not  catastroph ical  ly  large,  the  real  part  of  the  index  of  refraction 
(Equation  2-184)  can  be  used  to  write  a  simple  expression  for  e0: 
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where  x  is  the  RF  wavelength,  R  is  the  distance  along  the  LOS  from 
transmitter  to  receiver,  and  rg  is  the  classical  electron  radius 
( 2  . 8 1 84 E - 1 3  cm) . 


A  similar  situation  exists  for  the  time  delay.  The  equation  for 
the  group  path  length  is  the  same  as  that  for  the  phase  path  length  except 
that  the  term  involving  the  electron  content  is  reversed  in  sign.  This 
leads  to  the  following  equation  for  time  delay  along  the  refracted  paths: 
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where  trf  is  the  time  delay  along  the  LOS  (including  electron  content 
thereon),  and  g^  is  given  by  a  summation  identical  to  that  of  Sj  (Equation 
2-194)  except  that  the  second  term  (involving  the  electron  content)  is 
reversed  in  sign.  Similarly  to  9  ,  td  may  be  written  in  the  following 
simpl i f ied  form: 
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2.5.4  Model  Implementation  Comments. 

Implementation  of  the  refraction  models  in  PRPSIM  follows  the 
descriptions  found  in  sections  2.5.1,  2.5.2,  and  2.5.3.  There  are,  how¬ 
ever,  a  few  details  and  assumptions  worth  mentioning. 

Typical  applications  of  the  multiscreen  modeling  require 
environment  data  at  a  set  of  points  along  the  LOS  only  within  a  spatial 
region  where  the  refractive  properties  are  not  negligible.  The  physical 
path  end  points  may  be  far  removed  from  the  intervening  refractive 
region.  The  multiscreen  method,  by  virtue  of  the  boundary  condition  that 
the  propagation  paths  must  connect  the  physical  path  end  points,  in  effect 
averages  the  properties  of  the  first  and  last  phase  screens  over  the  path 
segments  connecting  the  associated  end  points.  Depending  on  the  distances 
involved  and  the  refractive  properties  of  the  first  and  last  screens,  it 
is  sometimes  appropriate  to  insert  an  additional  "null  screen"  (i.e.,  one 
which  has  zero  spatial  derivatives)  between  one  or  both  of  the  physical 
end  points  and  the  associated  closest  calculation  point.  This  procedure 
(implemented  in  subroutine  SCREEN)  assures  that  regions  of  negligible 
refraction  are  treated  as  such.  Subroutine  REFRAC  calls  subroutine  SCREEN 
once  for  each  propagation  path  to  initialize  the  phase  screens. 


The  multiscreen  method  implemented  in  PRPSIM  requires  the  values 
of  the  spatial  derivatives  of  Ng  (mean  electron  density)  in  the  two 
orthogonal  directions  (vertical  and  traverse  planes),  up  to  fourth  order, 
at  each  screen.  DATAPT  (see  section  3.4),  calculates  the  3  components  of 
the  electron  density  gradient  and  the  6  components  of  the  second  spatial 
derivative  of  the  electron  density  at  each  screen  using  first  and  second 
order  differencing  respectively.  However,  this  technique  is  inappropriate 
for  the  higher  order  (third  and  fourth)  spatial  derivatives  needed  by  the 
multiscreen  method.  The  following  two  techniques,  dependent  upon  the 
magnitudes  and  signs  of  the  mean  and  the  first  and  second  spatial 
derivatives,  are  used  to  estimate  the  third  and  fourth  spatial  deriva¬ 
tives.  For  notational  convenience  we  define  the  following  values: 
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If  f?  <_  f2/f  then  a  Gaussian  profile  aproximation  is  used  for  the  higher 
order  derivatives: 
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The  above  inequality  provides  that  the  variance  of  the  Gaussian  profile 
will  be  sensible.  If  this  inequality  fails,  the  power-law  profile 
approx imat  ion  is  used: 
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The  use  of  this  profile  is  necessary  when  the  known  second  derivative,  f2, 
is  so  large  in  the  positive  direction  that  neither  the  Gaussian  nor  the 
exponential  profiles  is  a  reasonable  form  to  assume. 

The  higher  order  mixed  partials  are  calculated  using  the  usual 
product  function  approximations : 
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The  foe  ising/defocusing  parameter  e  (see  Equations  2-201  through 
2-204)  is  chosen  to  be  consistent  with  a  1  m rad  displacement  of  the  LOS. 
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In  some  cases  near  the  "caustic"  reqion  (where  there  is  no  KF  propagation 
due  to  the  extreme  ionization),  the  path  to  the  receiver  may  propaqate  but 
the  one  displaced  in  the  retractive  qradient  direction  will  not  (because 
it  is  beyond  the  caustic).  In  such  cases  we  rather  arbitrarily  set  the 
one-dimensional  one-way  defocusing  ratio  in  this  direction  equal  to  a 
small  value  (0.1)  correspondi nq  to  appreciable  focusing  which  is  typical 
of  this  tyoe  of  situation.  Conversely,  it  is  possible  in  some  cases  to 
obtain  only  one  real  root  for  the  cubic  perturbation  equation  (Equation 
2-19?)  on  the  path  to  the  receiver,  but  three  real  roots  on  one  of  the 
displaced  paths.  In  these  cases  two  of  the  three  roots  do  not  correspond 
to  physically  reasonable  paths,  and  hence  there  is  no  problem  as  long  as 
one  is  careful  to  associate  real  main  paths  with  real  displaced  paths  when 
computing  the  focusing  effects. 
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SECTION  3 

INTERFACE  WITH  THE  PROPAGATION  ENVIRONMENT 

3.1  INTRODUCTION. 

The  SCENARIO  code  was  developed  to  provide  reliable  estimates  of 
the  phenomenoloqy  of  hiqh  altitude  nuclear  bursts.  Preexisting  research 
pbenomeno loqy  codes  were  capable  of  providing  the  necessary  data  for 
individual  situations,  but  they  are  much  too  expensive  and  difficult  to 
run  to  be  of  practical  use  in  analyzing  a  wide  variety  of  nuclear  scen¬ 
arios.  SCENARIO  is  thus  a  compromise  in  that  it  calculates  only  those 
environmental  parameters  which  have  the  most  significant  impact  on  the 
propagation  of  satellite  signals.  Although  the  SCENARIO  geometry  cannot 
provide  as  high  a  degree  of  resolution  at  early  times  as  the  research 
codes  provide,  it  is  adequate  to  model  a  very  broad  range  of  realistic 
nuclear  scenarios.  The  fundamental  objective  of  the  SCENARIO  code  is  to 
predict  the  evolution  of  the  power  spectral  density  function  of  the 
electron  density  fluctuations  over  a  CONUS  sized  region  on  a  time  scale  of 
hours  during  and  following  a  series  of  high  altitude  nuclear  events.  The 
purpose  of  this  section  is  to  provide  the  user  with  a  very  basic  descrip¬ 
tion  of  the  data  contained  in  the  files  generated  by  the  SCENARIO  code  and 
how  the  PPPSIM  code  inr.erfaces  with  them.  For  a  more  complete  description 
of  the  SCENARIO  code,  the  user  should  refer  to  the  reports  MRC-R-40412, 
MRC-R-S39n,  ar  MRC-R-994  14.  For  a  description  of  the  format  of  SCENARIO 
data  files,  refer  to  MRC-N-718R 1 s . 
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3.2 


RETRIEVAL  AND  STORAGE  OF  SCENARIO  DATA. 


The  storage  and  retrieval  of  MHD  environmental  data  read  from 
SCENARIO  data  files  requires  more  CPU  time  than  all  other  PRPSIM  opera¬ 
tions.  For  that  reason,  simulation  time  is  designated  to  be  the  outermost 
loop  in  the  hierarchy  of  PRPSIM' s  structure.  With  such  a  structure,  data 
for  a  single  grid  dump  time  need  be  stored  only  once  per  simulation, 
reqardless  of  the  number  of  links,  frequencies,  etc.  requested  for  that 
simulation  time.  Environmental  data  are  read  from  SCENARIO  files  and 
stored  in  memory  by  the  PROPEN  routine.  They  may  then  be  retrieved  from 
memory  by  the  DATAPT  routine.  PROPEN  is  called  once  by  PRPSIM  for  each 
simulation  time  specified  -  provided  that  the  simulation  time  does  not 
correspond  to  the  grid  dump  time  of  a  file  already  in  memory,  or  that  it 
is  not  bracketed  by  the  grid  dump  times  of  two  files  already  in  memory. 

In  order  for  PROPEN  to  operate,  it  must  share  the  following  data 
with  PRPSIM  in  common  storage:  (1)  the  number  of  SCENARIO  data  files 
available  to  the  simulation,  (2)  the  array  containing  the  names  of  those 
riles  (here  referred  to  as  the  "file  name  array"),  (3)  the  array  of  the 
initial  grid  times  in  those  files  (here  referred  to  as  the  "initial  time 

array"),  (4)  the  total  number  grid  times  in  the  designated  files  (a  physi¬ 
cal  file  may  contain  more  than  one  grid  dump),  (5)  the  array  of  all  grid 

times  in  the  designated  files  (here  referred  to  as  the  "grid  time  array"), 

(61  the  current  simulation  time,  and  (7)  the  grid  dump  time(s)  of  the  data 
currently  in  extended  memory.  (Note  that  the  number  of  times  in  the 
initial  time  array  eauals  the  number  of  names  in  the  file  name  array.) 

It  is  assumed  that  the  data  file  names  are  ordered  in  increasing 
initial  grid  time  and  that  the  data  within  individual  files  (should  the 
file  contain  morp  than  one  grid  time)  are  ordered  in  increasing  grid 
time.  Also,  the  initial  grid  time  of  a  file  must  be  larger  that  the  last 
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grid  time  of  the  previous  file.  It  is  the  resoons i hi  1  i ty  of  the  INDATA 
module  to  insure  that  these  reouirements  are  satisfied.  Upon  being 
called,  PROPEN  first  checks  to  determine  whether  the  current  simulation 
time  is  within  the  ranqe  of  qrid  dump  times  in  the  designated  files.  If 
it  is  not,  an  error  messaqe  is  printed  and  the  iob  is  aborted  immediately. 

All  file  manipulation  and  bookkeeping  are  accomplished  by  two 
pointers.  Upon  beinq  called,  PROPEN  sets  the  "file  pointer"  to  the  larqest 
time  in  the  initial  time  array  which  does  not  exceed  the  current  simula¬ 
tion  time.  If  the  file  pointer  chanqes  from  its  previous  value,  the 
currently  open  data  file  is  closed.  Similarly,  the  "qrid  pointer"  is  set 
to  the  larqest  time  in  the  qrid  time  array  which  does  not  exceed  the 
current  simulation  time.  A  file  inauiry  is  conducted  to  determine  if  the 
file  in  the  file  name  array  indicated  by  the  file  pointer  can  be  located 
on  disk.  If  not,  an  error  message  is  printed  and  the  job  is  aborted. 
Note  that  PROPEN  always  returns  to  this  point  in  the  code  whenever  the 
file  pointer  is  incremented.  If  the  file  indicated  by  the  file  pointer  is 
not  open,  it  is  opened  and  rewound.  Before  proceeding  further  with  this 
description  of  the  PROPEN  routine,  the  user  should  refer  to  MRC-N-718R15 
to  familiarize  himself  with  the  format  of  the  SCENARIO  data  files. 

A  SCFNARTO  loqical  file  consists  of  the  following  records:  (1) 
the  header  record,  (?)  the  burst  history  record,  (3)  the  qrid  qeometry 
record,  (A)  the  ambient  atmosphere  ionosphere  record,  (5)  the  active/ 
inactive  tube  mao  record,  (6)  the  first  plasma  tube  record,  (7)  the  second 
plasma  tube  record,  and  so  on.  Two  optional  records  may  follow  the  data 
records.  Neither  of  these  optional  records  is  used  by  PRPSTM,  and  they 
will  not  he  discussed  further  here. 

PPOPFN  first  attempts  to  mad  a  header  record  from  the  designa¬ 
ted  CC^JARIQ  file.  Should  this  be  an  end-of-file  (EOF)  marker,  PROPEN 
thm  attempts  to  read  the  next  record.  If  more  than  two  EHFs  are  read  in 
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succession,  PROPEN  assumes  that  it  has  reached  the  end  of  the  physical 
file.  It  then  closes  that  file,  increments  the  file  pointer  by  one,  and 
returns  to  the  point  indicated  above  to  inquire  about  the  existence  of  the 
next  file  in  the  file  name  array.  If  the  record  is  not  an  EOF,  PROPEN 
assumes  that  the  record  is  the  header  record  of  the  next  logical  file  and 
proceeds  to  read  that  record  and  to  store  the  information  contained  on 
it.  The  header  record  contains  all  the  information  that  PROPEN  needs  in 
order  to  be  able  to  interpret  the  information  contained  in  all  the  remain¬ 
ing  records  in  the  logical  file.  The  first  information  from  the  header 
record  to  be  used  by  PROPEN  is  the  qrid  dump  time  which  determines  whether 
the  data  from  this  qrid  dump  (logical  file)  will  be  stored.  To  determine 
whether  the  data  in  the  current  loqical  file  should  be  stored,  PROPEN 
compares  the  current  simulation  time  to  the  current  grid  dump  time  and  the 
time  in  the  grid  time  array  to  which  the  grid  pointer  points.  If  the 
simulation  time  equals  the  qrid  dump  time,  or  if  it  is  less  than  either 
the  indexed  array  time  or  the  following  array  time,  the  data  will  be  read 
and  stored  by  PROPEN;  otherwise,  PROPEN  merely  reads  over  all  the  records 
in  the  current  qrid  dump  with  dummy  read  statements  and  then  returns  to 
the  point  in  the  code  at  which  to  read  the  header  record  of  the  next  file 
(or  the  E0F[sj  at  the  end  of  the  file.)  Conversely,  if  the  data  were 
stored ,  and  if  the  grid  dump  time  equals  or  exceeds  the  simulation  time 
(so  that  no  additional  data  are  required  for  time  interpolation),  PROPEN 
returns  control  to  PRPSIM  after  readinq  the  last  record  in  the  file.  Note 
that  the  current  input  data  file  is  always  left  open  when  PROPEN  returns 
control  to  PRPSIM. 


The  contents  of  the  second  record,  the  event  history  record,  are 
read  and  printed  by  PROPEN  it  the  diagnostic  level  is  set  to  1  or  higher. 
Otherwise,  they  are  not  used  by  PRPSIM. 
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The  third  record  is  the  qrid  geometry  record  that  contains  all 
the  information  needed  by  the  path  intersection  routine  (PATHPT)  to  deter¬ 
mine  the  intersection  of  the  propagation  path  with  the  grid.  Because  the 
size  and  location  may  or  may  not  change  from  one  grid  dump  to  another, 
PROPEN  maintains  a  logical  array  (GRDFLG)  that  is  shared  with  PATHPT  to 
inform  it  of  any  changes  in  the  grid  geometry.  GRDFLG(l)  is  set  true  if 
the  grid  geometry  data  in  register  one  are  replaced  by  new  data,  and 
GRDFLG(2)  is  set  true  if  those  in  register  two  are  replaced.  They  are 
reset  to  false  by  PATHPT  when  it  processes  the  new  geometry  data. 
GRDFLG(3)  is  set  false  by  PROPEN  if  the  geometries  of  the  grids  for  the 
data  in  both  registers  are  identical;  otherwise  it  is  set  true. 

The  fourth  record  is  the  ambient  atmosphere  ionosphere  record 
which  contains  a  simple  (radial)  ambient  atmosphere/ionosphere  model  which 
the  data  retrieval  module  (OATAPT)  can  use  to  generate  data  for  points 
found  in  inactive  tubes.  Because  the  ambient  model  may  also  change  with 
time,  PROPEN  maintains  two  registers  in  which  these  data  are  stored. 

The  fifth  record  is  the  act i ve/i nact i ve  plasma  tube  mao  record. 
This  record  consists  of  a  sequential  list  of  indices  for  every  data  tube 
in  the  grid.  If  the  tube  is  active,  its  index  s  its  position  in  the 
sequence  of  data  tube  records;  if  it  is  inactive,  its  index  is  zero. 

The  remaining  records  are  the  (active)  plasma  data  tube 
records.  Current.lv,  files  are  written  with  one  tube  oer  data  record. 
PR OPEN  can  read  multiple  tube  records,  however,  provided  that  the 
necessary  information  is  included  on  the  header  record. 

Normally,  PROPEN  maintains  two  registers  in  which  to  store  the 
SCENARIO  data  so  that  data  can  be  interpolated  between  simulation  times 
which  do  rot  coincide  with  a  grid  dump  time.  Should  the  storaqe 
'•egu  i  rement  s  of  two  registers  exceed  the  storaqe  ranacity  available  to  you 


79 


I 


on  your  machine,  the  second  register  may  be  dispensed  with  --  provided 
that  no  time  interpolation  is  required.  This  is  accomplished  readily  by 
setting  the  value  of  the  parameter  ISTACK  to  1,  rather  than  the  usual 
value  of  2.  When  both  registers  are  used,  storage  of  the  SCENARIO  data 

alternates  between  them.  The  value  of  of  the  register  pointer  I  REG  (1  or 
2)  indicates  which  register  was  last  filled.  Each  register  has  its  own 
time  reqister  which  is  equated  to  the  current  grid  dump  time  when  data  are 
stored  in  the  register. 

On  virtual  storaae  machines  such  as  the  VAX/VMS  system,  data 

from  the  active  plasma  tubes  are  stored  in  a  large  array  which  is  accessed 
only  by  the  storage  routine  PROPEN  and  the  retrieval  routine  DATAPT.  On 
other  machines  such  as  the  CDC/Cyber  series,  the  data  are  stored  in 
extended  core  memory.  All  other  data  read  by  PROPEN  are  stored  in  arrays. 

Table  4.  SCENARIO  cell  quantities  stored  in  mass  storage  by  PROPEN 

Niodu  1  e . 


PARAMETERS 

UNITS 

I'D 

mean  electron  density  in  the  cell 

[cm-3] 

|  ( ?) 

standard  deviation  of  electron  density  in  the  cell 

[cm- 3] 

|  (  3) 

partial  time  derivative  of  mean  electron  density 

[cm-3  sec*1] 

(4) 

mean  electron  temperature  within  the  cell 

no 

(S) 

a  component  of  the  mean  ion  velocity 

[cm/sec] 

|(6) 

0  component  of  the  mean  ion  velocity 

[cm/ sec] 

|m 

*  component  of  the  mean  ion  velocity 

[cm/sec] 

|(8) 

mean  neutral  mass  density 

[gm/cm] 

|(<0 

_ 

transverse  outpr  striation  scale  size 

[cm] 
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3.3 


PROPAGATION  PATH  INTERSECTION  WITH  THE  SCENARIO  GRID. 


PATHPT  is  the  PRPSIM  routine  which  generates  the  set  of  analysis 
points  along  the  line  of  sight  (LOS)  connecting  each  pair  of  linked  ob¬ 
jects.  PATHPT  is  called  once  by  RFPROP  for  each  link  at  each  simulation 
time.  The  qeographic  coordinates  of  the  transmitter  (origin)  and  receiver 
(end)  are  passed  to  PATHPT  by  RFPROP  as  arguments.  PATHPT  then  proceeds 
to  calculate  the  number  of  cells  in  the  MHD  environmental  qrid  which  are 
intersected  by  the  LOS,  the  lengths  of  the  path  segments  that  are  defined 
by  the  intersections  of  the  LOS  and  the  cell  boundaries,  and  the  distance 
froin  the  LOS  origin  to  the  midpoint  of  each  of  the  segments.  The 
midpoints  of  the  segments  define  the  set  of  analysis  points  at  which 

RFPROP  calculates  the  integrands  needed  to  compute  the  various  signal 
parameters . 

A  logical  array  (GRDFLG)  is  shared  between  PATHPT  and  the  envi¬ 
ronmental  data  storaqe  routine  PROPEN  (described  in  section  3.2)  in  order 
to  avoid  unnecessary  work  by  PATHPT  when  the  grid  geometry  does  not  change 
between  successive  qrid  dumps.  In  qeneral ,  PROPEN  maintains  two  registers 
in  which  qrid  dumps  at  two  different  times  are  stored.  The  first  element 
of  the  loqical  array  is  set  true  by  PROPEN  if  (and  only  if)  the  qrid 
geometry  of  the  new  dump  to  be  stored  in  reqi^ter  one  differs  from  that 
nf  the  old  dump  (which  is  to  be  replaced)  in  reqister  one.  The  second 

element  is  analogous  and  refers  to  the  data  stored  in  the  second  reqis¬ 
ter.  If  either  of  these  two  loqical  flags  is  true,  PATHPT  calculates  the 
coordinates  of  the  qrid  boundaries  from  the  cel!  widths  and  the  cell 
center  coordinates  provided  by  PROPEN  and  sets  the  flaq  false.  Otherwise, 
it  calls  the  subroutine  XIT3X2  which  is  responsible  for  the  calculation  of 
tne  set  of  tne  path  segment,  lenqths  and  the  distance  to  the  midpoint  of 
tne  first  segment.  The  third  element  of  the  logical  array  is  set  true  by 
/•’:  N  if  tne  grid  geometries  of  the  dumps  in  the  two  registers  differ 

from  one  another.  If  this  flag  is  set,  and  if  time  interpolation  is 

required,  xlrjX2  is  called  twice,  otherwise  it  is  called  only  once.  (The 
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manner  in  which  the  two  sets  are  merged  to  form  a  single  set  of  inter¬ 
sections  will  be  described  after  we  describe  the  operation  of  the  sub¬ 
routine  X1T0X2.) 

The  first  task  of  X1T0X2  is  to  determine  the  maximum  and  minimum 
coordinates  of  the  LOS  along  each  of  the  three  grid  directions,  a,  8,  and 
<t> .  If  these  lie  outside  the  range  of  the  grid  boundaries  along  any  of  the 
three  grid  directions,  a  null  LOS-grid  intersection  is  obvious,  and  X1T0X2 
returns  control  to  PATHPT  with  the  number  of  intersections  set  to  zero. 
If  the  respective  maxima  and  minima  lie  within  the  limits  of  the  grid 
along  each  of  the  three  grid  directions,  X1T0X2  proceeds  to  calculate  the 
distance  to  the  intersection  of  the  LOS  with  each  of  the  boundary  surfaces 
of  the  grid.  Note  that  when  we  calculate  the  distances,  we  may  deal  with 
the  intersections  with  each  of  the  three  sets  of  surfaces  (i.e.,  those 
surfaces  with  constant  a,  8,  and  <t> ,  respectively)  independently  of  the 
other  two.  After  computing  all  the  intersections,  it  is  only  necessary  to 
combine  the  three  sets  of  distances  so  that  they  are  arranged  in  order  of 
increasing  distance  from  the  transmitter. 

The  set  of  surfaces  of  constant  are  easiest  to  visualize: 
tney  are  simply  a  set  of  planes  with  a  common  line  of  intersection  which 
coincides  with  the  dipole  axis.  The  8  surfaces  are  curves  with  an  axis  of 
rotational  symmetry  which  coincides  with  the  dipole  axis,  and  the  a  sur¬ 
faces  are  curves  th  a  plane  of  symmetry  which  coincides  with  the  dipole 
equator  (see  figures  5  and  6).  Let  us  define  a  Cartesian  set  of 
coordinates  with  the  z-axis  aligned  along  the  dipole  axis  (pointinq  toward 
the  north),  and  the  v-axis  defined  by  the  cross  product  of  the  rotation 
axis  and  the  dipole  axis.  The  components  of  any  point  alonq  the  LOS  may 
be  computed  as 
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gure  5.  Cross-sectional  view  of  a  phi-plane  in  the  dipole  field  aligned 
coordinate  system  of  a  uniform  SCENARIO  grid.  Grid  cell 
boundaries  are  formed  by  surfces  of  constant  beta  (solid  lines) 
and  constant  alpha  (dashed  lines). 


Figure  6.  Cross-sectional  view  of  a  phi-plane  in  the  dipole  field  aligned 
coordinate  system  of  a  “staggered"  SCENARIO  grid.  Grid  cell 
boundaries  are  formed  by  surfaces  of  constant  beta  (solid  lines) 
and  constant  alpha  (dashed  line  segments). 


x  =  x0  +  at  (3-1) 

y  -  y0  +  bt  (3-2) 

z  =  z0  +  ct  (3-3) 

where  t  u  the  distance  from  the  origin  to  the  point  along  the  LOS;  x0, 
y0,  and  z0  are  the  coordinates  of  the  origin;  and  a,  b,  and  c  are  the 
direction  cosines  of  the  LOS  along  those  axes,  respectively, 

a  =  (VV^LOS  (3-4> 

b  =  (VV^LOS 

c  =  (ze-zo)7tLOS  (3‘6) 

where  xg ,  yg ,  and  zg  are  the  coordinates  of  the  end  point  of  the  LOS,  and 

lL0S  =  [<xe-x„)2  +<V*o>2  +  <ze-zo>2l‘/2  (3-7) 

is  the  length  of  the  LOS. 

To  compute  the  distances  to  the  intersections  with  the  set  of  * 
planes  with  longitude  <t>^ ,  we  use  the  relation 

,  y* 

*•  =  tan-1  f— -)  (3-8) 

1  x . 

l 

where  i  ranges  between  1  and  the  number  of  <f>  planes  defined  in  the  grid. 
Then,  from  Eguations  3-1,  3-2  and  3-8  we  find  the  oistances  to  the 
inter  sect i ons 

y  cos^i  -  x  sin*. 

t;  =  - — - - - -  (3-9) 

1  asm*.  -  bcos*. 
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Note  that  if  the  LOS  (when  extended  to  infinity  in  both  directions)  inter¬ 
sects  the  dipole  axis,  there  are  no  intersect  ions  of  the  LOS  and  the  <t> 
pi anes . 


Equation  3-9  yields  the  intersection  of  the  LOS  with  every  <j> 
plane  in  the  grid.  A  real  grid  is  bounded  in  a  and  8  (as  well  as  if)  so 
chat  we  must  test  each  of  the  points  to  determine  whether  it  lies  within 
these  limits.  Calculate  the  alpha  coordinate  of  the  ith  point  using 


cos  a  ■ 


z 


2 

e 


(3-10) 


and  the  Deta  coordinate  using 


s  i  n 2  B  . 


Re 

(x  2  +  y2  +  z2) 3/2 


(3-11) 


where  Rg  is  the  radius  of  the  earth,  and  x-,  y.,  and  z-  are  obtained  by 
inserting  t-  from  Equation  3-9  into  Equations  3-1,  3-2,  and  3-3.  If 
either  of  these  values  falls  outside  the  boundaries  of  the  grid,  the  ith 
intersection  is  excluded  from  the  set  of  grid  intersections. 

Equations  3-10  and  3-11  may  be  used  to  calculate  the  inter¬ 
sections  of  the  LOS  with  the  a  and  P  surfaces,  respectively.  The 
distances,  t.  ,  to  the  intersections  with  the  R  surfaces  are  given  by  all 
the  real  roots  of  the  sixth  order  polynomial 

r(x,  +  at.)9  +  fyo+bt.)2  +  (z0+ctj)?l3  =  R|r(xQ+atj)2  +  (y0+bt  ^  )  2 1  2/s  in  48j. 

(3-12) 

such  that  0  '  t.  •  tLQ<-  where  t^  is  the  length  of  the  LOS,  and  j  ranges 
over  the  number  of  r  surfaces  in  the  grid.  Similarly,  the  distances  to 
the  intersect  ions  with  the  n  surfaces  are  given  by  the  real  roots  of  the 
sixth  order  polynomial 
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r(x0+atk)?  +  (y0+btk)?  +  (z0+Ctk)213  =  Re4(Vctk)2/cos\  (3-13) 


such  that  0  <  tk  <  t^^  where  t^^  is  the  length  of  the  LOS,  and  k  ranges 
over  the  number  of  n  surfaces  in  the  grid.  The  SCENARIO  grid  is  compli¬ 
cated  by  the  fact  that  the  a  coordinates  are  "staggered"  along  different  R 
coordinates.  Specifically,  for  smaller  values  of  R,  the  n.  coordinates 
decrease  toward  smaller  values,  but  the  number  of  a  coordinates  is  the 
same  for  each  of  the  B. 

The  number  of  real  roots  of  Equation  3-12  on  the  interval 
[0,tL0S  ]  is  determined  by  applying  Sturm's  theorm.  The  first  of  these 
roots  is  then  calculated  precisely  by  applying  Laguerre's  method  for  root 
extraction  with  a  "rough"  initial  estimate  on  the  interval  [0,tLg^].  (For 
those  rare  instances  in  which  Laguerre's  method  diverges,  Sturm's  method 
can  be  used  to  obtain  a  better  initial  estimate  of  the  root,  so  that 
Laguerre's  method  can  be  applied  again  to  converge  to  the  root.)  If  more 
than  one  root  (intersection)  exists,  the  polynomial  is  "deflated"  by 
dividing  out  the  previously  calculated  root,  and  the  root  of  the  deflated 
polynomial  is  calculated  as  before.  This  process  is  continued  until  all 
the  (up  to  four)  real  roots  on  the  interval  [0 , L^qs 1  are  ca1culated.  Each 
of  these  distances  must  be  checked  by  calculating  (x . ,  y . ,  z.)  for  each  t. 
(using  Equations  3-1,  3-2,  and  3-3)  and  inserting  these  values  into  Equa¬ 
tions  3-8  and  3-10  to  insure  that  the  corresponding  <*  and  .f>  coordinates 
lie  within  the  boundaries  of  the  grid. 

The  roots  of  Equation  3-13  are  calculated  in  analogous  manner.  Because 
the  grid  is  staggered,  however,  it  is  convenient  to  analyze  individually, 
each  seqment  of  the  LOS  between  the  intersections  with  adjacent  a  sur¬ 
faces.  If  one  (or  both)  end(s)  of  the  LOS  lie(s)  within  the  grid  bound¬ 
aries,  it  is  necessary  to  add  one  (or  two)  segments  between  the  beginning 
and/or  end  of  the  LOS  and  the  first  and/or  last  intersected  R  surface. 
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The  distances  to  the  intersections  between  each  of  these  segments  and  the 
i  surfaces  in  the  "tube"  in  which  the  segment  lies  are  calculated  by  find¬ 
ing  the  real  roots  of  Equation  3-13  on  the  interval  [0,t-]  where  t-  is  the 

J  J 

length  of  the  jth  segment.  (The  Xq,  y^,  and  Zq  of  Equation  3-13  must  be 
replaced  by  the  Cartesian  coordinates  of  the  beginning  of  the  line  seg¬ 
ment.)  The  roots  of  this  sixth  order  polynomial  are  computed  in  the  same 
manner  as  those  of  Equation  3-12. 

When  the  three  sets  of  distances  to  the  intersections  of  the  LOS 
with  each  of  the  three  sets  of  surfaces  have  been  computed,  X 1 T 0X2  com¬ 
bines  all  three  sets  and  reorders  them  according  to  increasing  distance 
from  the  origin.  The  parameters  that  are  passed  back  to  PATHPT  are  the 
differences  between  distances  to  adjacent  intersections  (i.e.,  the  set  of 
path  segment  lengths,  half  the  sum  of  the  distances  to  the  first  two 
intersections  (i.e.,  the  distance  to  the  center  of  the  first  segment),  and 
the  number  of  intersections  minus  one  (i.e.,  the  number  of  segments). 

Because  the  grid  is  bounded  by  curved  surfaces  in  a  and  B,  it  is 

possible  tnat  the  LOS  may  exit  the  grid  across  a  concave  surface  and  then 

reenter.  In  such  cases,  X1T0X2  generates  a  path  segment  length  for  the 
extra-grid  portion  of  the  LOS  as  though  it  were  enclosed  by  an  ordinary 
grid  cell.  The  data  retrieval  routine  DATAPT  recognizes  this  situation  and 
returns  null  data  (zeros)  for  such  a  segment.  Path  segments  outside  the 

grid  but  not  bracketed  by  segments  inside  the  grid  are  not  generated  by 

X1T0X2 . 


When  time  interpo 1  at  ion  between  two  grids  with  different  geom¬ 
etries  is  required,  X1T0X2  must  be  called  by  PATHPT  twice  to  generate  two 
different  sets  of  segments  along  the  path.  These  must  be  merged  into  a 
single  set  of  segments  and  distances  to  be  used  by  RFPROP.  The  simplest 
means  of  doing  this  is  just  to  combine  both  sets  of  inter  sect  ions ,  and 
then  determine  the  set  of  segment  lengths  between  these  intersections.  If 
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the  number  of  segments  through  the  two  grids  is  Nl  and  N2,  then  the  number 
of  combined  segments  is  just  N1+N2+1.  Some  of  these  may  not  be  signifi¬ 
cant  segments,  however.  Consider  the  case  where  the  grids  at  the  two 
times  are  identical  but  the  second  is  shifted  a  minute  linear  distance 
relative  to  the  first.  (By  minute,  we  mean  that  the  distance  of  the  shift 
is  small  when  compared  to  any  of  the  segment  lengths  along  the  LOS  through 
either  of  the  grids.)  As  one  can  easily  visualize,  the  effect  of  combin¬ 
ing  the  segments  at  the  two  times  in  the  manner  described  above  would  be 
to  insert  minute  segments  of  length  equal  to  the  shift  between  the  origi¬ 
nal  segments.  In  such  a  case,  RFPROP  would  have  to  perform  calculations 
at  twice  the  number  of  points  it  actually  needs.  Our  method  of  eliminat¬ 
ing  what  we  consider  to  be  insignificant  segments  from  the  combined  set  is 
to  merqe  a  segment  with  its  shortest  neighbor  if  its  length  is  less  than 
one  fourth  the  length  of  its  shortest  adjacent  neighbor.  Also  all  seg¬ 
ments  shorter  than  5  m  are  merged  with  their  shortest  neighbor. 

Finally,  before  returning  control  to  RFPROP,  PATHPT  calculates 
the  array  of  distances  from  the  origin  to  the  midpoint  of  each  of  the  path 
segments.  These  constitute  the  set  of  points  along  which  the  propagation 
integrals  will  be  calculated. 

3.4  CALCULATION  OF  PRIMITIVE  DATA  FOR  PROPAGATION  INTEGRALS. 

A  set  of  nine  SCENARIO  data  exists  for  every  active  cell  within 
the  plasma  grid.  Specifically,  each  of  these  data  represents  the  mean 
value  (and,  in  the  case  of  the  electron  density,  the  standard  deviation 
about  the  mean)  of  a  particular  variable  within  the  cell.  For  a  particu¬ 
lar  PRPSIM  simulation  of  the  propagation  of  a  signal  between  a  transmitter 
and  receiver,  a  number  of  points  along  the  propagation  path  is  generated 
by  the  routine  PATHPT.  In  theory,  this  set  of  points  is  the  optimum  set 
for  the  particular  plasma  grid  because  of  the  one-to-one  correspondence 
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between  the  number  of  points  and  the  number  of  grid  cells  which  the  propa¬ 
gation  path  intersects.  (When  time  interpo 1  at  ion  between  two  incongruent 
grids  is  required,  the  optimum  set  of  points  is  harder  to  define,  but  that 
problem  is  beyond  the  scope  of  this  section.)  The  basic  task  of  the  data 
retrieval  routine  ( DATAPT )  is  to  provide  the  integration  module  RFPROP 
with  a  set  of  data  for  each  of  these  points  along  the  path.  The  simplest 
means  of  doinq  this  is  to  equate  the  data  at  each  analysis  point  to  the 

mean  values  for  the  cell  in  which  it  is  located.  DATAPT  utilizes  a  method 
which  somewhat  more  accurately  represents  the  spatial  variations  of  the 

variables  in  the  grid  by  means  of  a  three  dimensional  interpolation  of  the 
data  from  the  coordinates  of  the  nearest  grid  cell  centers  to  the  coordi¬ 
nates  of  the  analysis  point. 

DATAPT  operates  in  either  of  three  computational  modes.  For  the 
basic  mode  in  which  the  spatial  derivatives  of  the  electron  density  are 
computed,  data  for  four  cells  in  each  of  twelve  tubes  surrounding  the 

point  are  fetched.  The  data  at  the  eight  cell  centers  adjacent  to  the 

point  are  used  for  spatial  interpolation,  and  the  electron  densities  from 
the  24  "next  nearest"  neighbors  are  used  to  calculate  the  density  gradient 
and  second  derivatives.  For  the  second  mode  in  which  the  derivatives  are 
not  computed,  data  from  the  eight  nearest  cell  centers  only  are  fetched 
and  i nterpo 1 ated .  The  third  mode  is  used  for  antenna  temperature  calcula¬ 
tions  and  is  similar  to  the  second  mode,  but  only  a  truncated  data  set  is 
returned  to  RFPROP.  In  order  to  minimize  the  number  of  data  fetches, 

DATAPT  maintains  an  index  of  the  subtubes  for  which  data  are  currently 

available  in  the  local  buffer.  (Here,  a  "subtube"  refers  to  4  [mode  1]  or 
2  [modes  2  and  3]  adjacent  cells  in  a  tube.)  For  subsequent  points  along 
the  path,  the  set  of  subtubes  for  which  data  are  needed  is  compared  to 
those  already  in  the  local  buffer.  The  number  of  subtubes  which  must  be 
fetched  equals  the  number  of  subtubes  in  the  buffer  which  are  no  longer 

needed,  so  that  the  former  may  be  written  over  the  latter,  and  the  index 

updated.  Figure  7  shows  a  schematic  representat ion  of  the  portion  of  the 
SCENARIO  grid,  relevant  to  a  DATAPT  operation  at  one  point. 
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Figure  7.  Cartesian  representation  of  the  portion  of  the  SCENARIO  grid 

used  for  a  OATAPT  calculation  at  a  single  point.  The  data  from 
the  8  nearest  cell  centers  (solid  circles)  are  used  to  inter¬ 
polate  data  to  the  coordinates  of  the  point.  The  remaining  24 
•next-nearest"  neighbors  (open  circles)  are  used  to  approximate 
the  mean  electron  density  gradient  and  the  spatial  second 
derivatives.  (Note  that  the  grid  shown  is  not  "staggered.") 


91 


Thp  data  storage  routine  PROPEN  insures  that  the  current  simula¬ 
tion  time  equals  tne  grid  time  of  one  of  the  grid  data  dumps  currently  in 
extended  memory  or  that  the  grid  times  of  the  stored  data  bracket  the 
simulation  time.  In  the  latter  case,  DATAPT  must  perform  its  calculations 
at  each  point  twice  using  the  data  at  both  grid  dump  times.  A  single  data 
set  is  then  calculated  from  the  two  sets  by  in  terpo  1  at  ion  in  time. 
(Obviously,  two  local  buffers  are  required  to  accomodate  the  data  from  all 
subtubes  at  both  grid  dump  times.) 

In  addition  to  interpolation  of  SCENARIO  data,  DATAPT  calculates 
several  of  its  own  parameters.  These  parameters  calculated  by  DATAPT  and 
interpolated  from  SCENARIO  data  are  summarized  in  Table  5. 


Table  5.  Summary  of  parameters  calculated  or  interpolated  by  DATAPT 
submodule . 


I  =.*  interpolated  from  SCENARIO  data 

C  =>  calculated  by  DATAPT 

C/I  =/■  calculated  using  the  interpolated  SCENARIO  data 

Mode  =  1/2  =>  ordinary  calculation  with/without  local 

electron  density  derivatives 

Mode  =3  antenna  temperature  calculation 

PARAMETER 

UNITS 

MODE 

mean  local  electron  density 

[cm-3] 

I 

1 

2 

0 

electron  temperature 

[  °K  ] 

I 

1 

2 

3 

electron  density  standard  deviation 

[cm-3] 

I 

1 

2 

3 

rms  electron  density 

[cm- 3] 

C/I 

1 

2 

3 

electron-neutral  collision  frequencies 

[sec-  !] 

C/I 

1 

2 

3 

electron-ion  collision  frequency 

[sec-1] 

C/I 

1 

2 

3 

3  components  of  the  ion  velocity 

[cm  sec- 

5] 

I 

1 

2 

partial  time  derivative  of  electron  density 

[cm- 3  sec 

-1] 

I 

1 

2 

?  transverse  outer  striation  scale  sizes 

[cm] 

I 

i 

2 

orientation  of  outer  transverse  scale  sizes 

[rad] 

C 

1 

2 

parallel  outer  striation  scale  size 

[cm] 

c 

1 

2 

isotropic  inner  striation  scale  size 

[cm] 

c 

1 

2 

3  components  of  the  magnetic  field 

[Gauss] 

c 

1 

2 

3  components  of  electron  density  gradient 

[cm-] 

C/I 

1 

6  components  of  second  derivative  of  density 

[cm- 

C/I 

1 
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The  arqument  list  for  DATAPT  consists  of  (1)  the  altitude  of  the 
point,  (2)  its  colatitude,  (3)  its  longitude,  (4)  the  number  of  propaga¬ 
tion  frequencies,  (5)  the  array  of  propagation  frequencies,  (6)  the 
sequential  number  of  the  point  on  the  path  (used  for  diagnostic  printout 
only),  and  (7)  a  logical  array  which  specifies  the  operational  mode  for 
the  path. 


The  qeomagnetic  coordinates  (geomagnetic  colatitude  and  longi¬ 
tude,  respectively)  of  the  point  are  calculated  from  its  geoqraphic 
coordinates  using 

=  arcco s  [coseDcos0Q  +  sineDsin90cos(4>D-4>0;  ]  (3-14) 


-sine0  sin(<t>D-<t>0) 

<t>  =  arctan  [ - 

m  cos9r)Sin0o  co$(<j>Q-ij>0 )  -  sin9QCOs8Q 


(3-15) 


where  4-j  is  the  aeograohic  colatitude  of  the  dipole  axis,  <j>p  is  the  geo¬ 
graphic  longitude  of  the  dipole  axis,  e0  is  the  qeographic  colatitude  of 
the  point,  and  t>0  is  the  geographic  lonqitude  of  the  point.  The  dipolar 
(a,  3,  t)  coordinates  are  then  calculated  from  the  geomagnetic  coordinates 
using 


COS  a 


s  i  n?, 


R 

R 

r_e 

'r 


2 

COS0 

m 

1/2 

sine 

m 


■t 


m 


(3-16) 


(3-17) 

(3-18) 


where  Rp  is  the  radius  of  the  earth  [cm],  and  R  is  the  radial  distance 
lot]  to  the  point  from  the  center  of  the  earth. 
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Given  the  dipolar  coordinates  of  the  point,  DATAPT  proceeds  to 
calculate  the  indices  (I^.I^I  )  of  the  cell  containing  the  point  by 

searching  through  the  arrays  of  grid  coordinates  read  from  the  grid  geo¬ 
metry  record  and  determining  the  pairs  of  grid  coordinates  which  bracket 
the  coordinates  of  the  point.  DATAPT  uses  the  convention  that  the  index 
alonq  a  given  axis  equals  the  sequential  number  of  the  cell  center  with 
the  largest  coordinate  which  does  not  exceed  the  coordinate  of  the  point 
on  that  axis.  The  geometry  of  the  SCENARIO  grid  is  complicated  by  the 
fact  that  the  t  coordinates  along  a  field  line  are  "staggered"  --  that  is, 
a  specific  o  index  corresponds  to  different  a  coordinates  on  different 
field  lines  (different  R  indices).  DATAPT  calculates  four  a  indices  for 
each  point  along  its  field  line  ( R  coordinate):  one  each  for  the  two  R 
indices  above  the  point  and  for  the  two  S  indices  below  the  point.  (Cur- 
reitly,  the  ot  index  of  a  point  in  the  SCENARIO  grid  ranges  between  n  and 
11,  and  the  Q  and  f  indices  range  between  0  and  3?.)  A  zero  index  indi¬ 
cates  that  the  point  lies  below  the  lowest  cell  center  on  that  axis,  and 
the  maximum  value  indicates  that  it  lies  above  the  highest  cell  center. 

Denote  the  indices  of  the  cell  containing  a  point  as  I I g,  and 
I  DATAPT  calculates  the  tube  indices  of  the  four  tubes  surrounding  the 
point  and  the  indices  of  the  eight  additional  "next-nearest  neighbor" 
tubes  (for  mode  1)  using 

"tube  =  *  Nx*(JS  -  »  <3-19> 

where  is  the  total  number  of  tubes  in  the  grid  along  the  R  axis  (read 
from  the  header  record).  ,10  ranqes  between  max(I  -1,1)  and  min(I^+2,Nx  ) 
for  mode  1,  and  between  max(I,,l)  and  min(Ic+l,Nx  )  for  modes  2  and  3; 

whereas  d  ranges  between  maxfl  -1,1)  and  min(I  +2,N  )  for  mode  1  and 

t?  <p  y 

between  max(I+,l)  and  min(I  +1,N  )  for  modes  2  and  3.  (The  min  and  max 

j  j  y 

functions  insure  that  DATAPT  does  not  attempt  to  retrieve  data  for  tubes 
which  do  not  exist.)  Note  that  the  four  "corner"  tubes  with  indices 
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( I  .,-1 , 1  -1 ) ,  (f.,-1 , 1 1+?) ,  (1 ,+2, 1  ^-1),  and  (l  ^2,1  s+2)  are  not  needed  to 
calculate  the  density  derivatives  and  are  not  stored.  A  schematic  diagram 
of  the  relevant  portion  of  the  grid  for  a  single  DATAPT  calculation  is 
depicted  in  Figure  7.  In  order  to  reduce  the  amount  of  data  transfer, 
DATAPT  does  not  fetch  the  contents  of  an  entire  grid  tube.  Since  we  need 
data  from  only  4  cells  for  mode  1  (at  2  cell  centers  above  the  point  and  2 
cell  centers  below  the  point)  or  2  cells  for  modes  2  and  3  (1  cell  center 
above  and  below  the  point),  DATAPT  fetches  data  in  units  of  these 
"subtubes."  (Because  of  the  staggering  of  the  grid,  it  is  important  to 
remember  that  the  a  indices  of  the  various  subtubes  will  differ  from  one 
another,  in  general.) 

DATAPT  then  checks  the  index  of  each  of  the  12  (or  4)  neighbor¬ 
ing  tubes  in  the  act i ve/ inact i ve  tube  map  index.  If  the  active/inactive 
tube  index  is  zero,  the  tube  is  inactive,  and  DATAPT  fetches  ambient  data 
for  the  the  locations  of  the  4  (or  2)  bracketing  cell  centers  and  stores 
them  in  the  appropriate  location  in  the  local  buffer.  If  the  active/ 
inactive  tube  index  for  the  tube  is  non-zero,  DATAPT  fetches  data  for  the 
bracketing  cell  centers  and  inserts  them  into  the  local  buffer.  From  this 
point  on,  DATAPT  makes  no  distinction  between  data  from  "active"  or 
"ambient"  culls. 


DATAPT  uses  the  four  subtubes  immediately  adjacent  to  the  point 
for  spatial  interpolation  of  SCENARIO  data  to  the  point.  For  each  of  the 
four  tubes,  we  interpolate  (log-linear)  each  datum  between  the  cell 
centers  with  u  coordinates  bracketing  that  of  the  point.  Then  we 
interpolate  in  0  to  yield  two  values  for  each  datum  at  the  two  f 
coordinates  bracketing  the  point,  and  finally,  we  interpolate  in  f>  to 
yi°ld  the  value  of  each  datum  interpolated  to  the  coordinates  of  the 
point.  In  order  to  compute  the  density  gradient,  we  calculate  and  store 
the  logarithm  of  the  electron  density  at  each  of  the  8  cells  surrounding 
the  point  and  at  the  24  'next-nearest '  cell  centers.  We  calculate  the 
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first,  and  second  derivatives  of  loq(Np)  in  a  alonq  each  of  the  four 

tubes  immediately  surroundinq  the  point.  (The  loqarithm  of  the  mean 

density  is  used  because  loq-linear  interpolation  of  the  density  more 
accurately  represents  the  denisty  variation  in  space  over  the  distances  in 
question.)  To  approximate  the  derivatives  alonq  each  tube,  let  f^ 
denote  the  values  of  the  loq  density  at  the  4  points  a-  alonq  the  a 

axis  with  a,  -  a,  <  a  ■'  a ,  <  a.,  where  a  is  the  a  coordinate  of 
the  point  at  which  the  derivative  is  to  be  calculated.  (See  Fiqure  7.) 
Owing  to  the  "staqqerinq"  of  the  a  coordinates,  the  a.  set  will  not, 
in  qeneral,  be  the  same  4  values  on  the  4  tubes.  Approximate  the  first 
derivatives  at  the  three  intermediate  points  halfway  between  a  3  and  a  2,  a2 
and  a3,  and  a 3  and  a4,  respectively: 


h±k 

a  2-a  j 


at  a21  =  ~  (a2+a  j) 


(3-20) 


,  =  f3~f2 

a3-a2 


at  a  32  =  —  ( a  3+a  2) 


(3-21 ) 


t,'t3 


at  a  43  =  _  (a4+a  3) 


( 3-22) 


Approximate  the  second  derivatives  at  the  2  points  near  a2  and 
a  3,  respectively: 


d  3  2~d  2  l 
a  3  2“ a  2 1 


at  a32i  -  —  (a32+a2i) 


(3-23) 


diQ'd  32 

a  43- a  32 


at  a432  -  —  (a  43+a  32) 


(3-24) 


The  value  of  the  first  derivative  at  the  point  a  =  a  is  then 

p 

approximated  by  linear  interpolation  between  the  first  differences  at  the 
two  bracketing  points: 


(3-25) 


3f 

3a 


d,,  +  j  i2-.1!?!  (a  -a  )  =  dd7,  +  dd  ,  (a  -a  ) 

21  a32-a21  0  21  21  321  D  2l; 


i  f  a  <  a 

D  -  32 


—  =  d32  +  K-3  )  =  dd32  +  dd4  32 

da  ^4  3*^32  ^  3  ^ 


i  i  a_  >  a 


32 


(3-26) 


In  addition  to  these  four  values,  the  derivatives  on  the  eiqht 
outer  tubes  are  calculated  using  the  two  values  bracketinq  the  a  coordi¬ 
nate  of  the  ooint.  These  values  are  stored  to  be  used  later  to  calculate 
the  mixed  second  derivatives  in  a  and  3.  and  in  a  and  a>. 


The  value  of  the  second  derivative  in  a  at  a  =  a  is  aporoxi- 

P 

mated  by  linear  interpci ation  between  the  second  differences: 

■lit  ,  dd„i  ♦  jl"u:dd.32i  (ap.<32i)  n-?7) 

)a  3432-332! 

To  calculate  the  derivatives  at  the  ooint  ( a  , 3  , t>  ),  we  linear- 

DPP 

ly  interpolate  each  pair  cf  values  in  3  to  yield  two  values  at  the  <j> 
coordinates  bracketinq  the  point.  Then  we  interpolate  in  $  between  these 
two  values  to  yield  the  final  approximation  to  the  derivatives  at  the 
coordinates  of  the  point. 


Calculation  of  the  derivatives  in  the  3  and  f  directions  is 
analogous  to  that  in  the  a  direction,  but  it  is  complicated  by  the  fact 
that  the  arid  is  "staggered."  Whereas  all  points  with  common  3  and  4> 
indices  have  the  same  3  and  coordinates,  the  same  is  not  true  for  the  a 
indices  and  coordinates.  A  set  of  points  with  the  same  a  and  3  indices, 
for  example,  do^s  not  constitute  a  "tube"  along  which  all  the  points  have 
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the  same  i  and  3  coordinates.  The  simplest  means  to  circumvent  this 
difficulty  is  to  interpolate  linearly  in  a  between  the  cell  centers  brack¬ 
eting  the  i  coordinate  of  the  point  along  each  (3,$)  tube  to  yield  a  grid 
of  values  for  the  function  on  the  " -x  surface"  containing  the  point.  The 
derivatives  in  3  can  be  calculated  each  of  the  two  lines  with  $  coordi¬ 
nates  bracketing  that  of  the  point.  The  two  sets  of  values  are  then 
interpolated  in  ,*>  t.o  yield  the  approximate  derivatives  in  3  at  the  noint. 
The  derivatives  in  *,  are  obtained  in  the  same  manner  t  /  i  nt  erpol  at  i  on  in 


The  mixed  second  derivatives  in  1  and  ?,  and  in  1  and  >  are 
calculated  in  exactly  the  same  manner  as  the  first  derivative  of  loq  (Np) 
using  the  derivatives  in  a  at  the  twelve  qrid  points  calculated  and  stored 
above.  The  mixed  second  derivative  in  B  and  *,  is  caluclated  using  the  six 
first  derivatives  in  j>,  plus  the  two  additional  differences  obtained  from 
the  other  four  tubes.  Differencinq  these  eight  values  in  3  yields  five 
estimates  of  3^f /  3 b  We  then  use  the  three  values  in  the  quadrant  of 
the  point  to  interpolate  for  the  value  at  the  point. 

If  the  simulation  times  does  not  equal  the  qrid  dump  time  for 
the  data  in  either  of  the  data  registers,  interpolation  in  time  of  all  the 
basic  qrid  data  is  required.  [n  this  case,  all  the  precedinq  PATAPT 
calculations  are  repeated  for  the  second  bracketing  grid  dump  time.  (It 
is  the  responsibility  of  the  routine  PRPPEM  to  insure  that  the  dump  times 
of  the  data  in  the  storaqe  registers  always  bracket  the  simulation  time.) 
Fach  pair  of  results  is  then  interpolated  (loq-1 inear)  to  the  simulation 
time. 


Thp  remainder  of  the  PflTAPT  code  is  devoted  to  the  calculation 
of  several  parameters  which  are  derived  as  functions  of  the  basic  data 
whose  calculation  is  described  above  and  as  functions  of  the  location  of 
the  point,  and  to  the  transformation  of  the  density  derivatives  from  the 
curvilinear  dipole  coordinate  system  to  a  ("artesian  svstpm. 
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The  mean  el ectron-neutral  momentum  transfer  collision  frequency 
[per  electron  per  second]  is  given  by 

ven  =  8. 14*  10!  2  on  Tj-64  [sec-1]  (3-28) 

where  k)  is  the  neutral  mass  density  [am  cm-1]  and  T.  is  the  electron 
temperature  [°K].  This  formula  is  derived  from  Equation  5-31  of  DNA- 
3R27T13  by  assuming  a  uniform  neutral  atmosphere  consisting  of  equal  parts 
of  monatomic  nitroqen  and  monatomic  oxygen  (that  is,  a  constant  mean 
molecular  weight  of  15  gm  mole-1). 

The  root  mean  square  [rms]  electron  density  is  given  by 

<N  >  =  /N2  +  a2  [cm-3]  (3-29) 

e _ _  e  N 

rms  e 

where  a is  the  electron  density  standard  deviation  [cm*3]  and  Ng  is  the 
mean  electron  density  [cm*3]. 

The  mean  electron-ion  momentum  transfer  collision  frequency  per 
electron  per  second  is  given  by 

T3 

v  .  =  1.8  <Np>  Tc1-5  In f 1 .25x10 16  —L)  [sec"1]  (3-30) 

'  rms  f2 

c 

where  ''N  '  is  the  rms  electron  density  [cm"3],  T.  is  the  local  electron 
e  rms  o  1 

temperature  [  K],  and  f  is  the  signal  carrier  frequency  [Hz].  (This 

equation  was  taken  from  page  5-1  of  DNA-3499H  1 7 . )  Note  that  an  array  of 

electron-ion  collision  frequencies  is  calculated  by  DATAPT  --  one  for  each 

carrier  frequency  in  the  simulation. 
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Collisions  between  ions  and  neutral  molecules  are  ignored  by 
DRPSI^  because  their  contribution  to  the  total  transfer  of  momentum  is 
negligible  in  comparision  with  collisions  between  electrons  and  neutral 
atoms  at  the  frequencies  and  altitudes  of  interest  here. 

DATAPT  sets  the  outer  striation  scale  size  parallel  to  the  field 
vector  to  a  constant  value  l_t  =  1.5E+7  cm.  This  value  insures  that 
the  ratio  of  the  parallel  scale  size  to  the  transverse  scale  is  in  the 
range  10  — IS. 

Currently,  the  SCENARIO  code  assumes  azimuthal  symmetry  about 
the  field  vector  and  generates  only  one  transverse  outer  striation  scale 
size.  DATAPT  assigns  an  equal  value  for  the  second  transverse  outer  scale 
size,  and  sets  the  orientation  of  the  axis  of  the  greater  scale  size  to 
zero . 


DATAPT  sets  the  isotropic  inner  striation  scale  size  to  a 
constant  value  l  =  1.0E+3  cm.  This  value  gives  a  minimum  outer  to  inner 
scale  size  ratio  of  approximately  1000. 

The  magnetic  field  vector  is  calculated  by  approximating  the 
earth's  field  by  that  of  a  geocentric  dipole  moment  of  8.1E+25  emu 
oriented  along  an  axis  through  78.6°  north  latitude  and  69.8°  west 

longitude  (epoch  1950;.  The  field  strength  is  calculated  as 

„  M  - 

J8|  =  —  /I  +  3cos2em  [Gauss]  (3-31) 

R 

where  is  the  dipole  moment,  R  is  the  radial  distance  to  the  point 

from  the  center  of  the  earth  [cm],  and  0m  is  the  geomagnetic  co  lati¬ 
tude  of  the  point.  By  definition,  the  dipole  field  vector  is  directed 

along  the  -  t  axis. 
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To  transform  the  components  of  the  ion  velocity,  magnetic  field, 
and  electron  density  gradient  from  the  field  aligned  coordinate  system 
(a,  R,  <t>)  to  the  local  tangent  plane  system,  calculate  the  inclination 
(dip  angle)  and  declination  of  the  magnetic  field  vector: 

PcosQ 

DIP  =  arctan  ( - _)  (3-32) 

sin  9^ 
m 


DEC  =  arctan  f 


sinn^  sine 
coso^  - 


o  sinQ0-»„) 
cosl0c°sflm 


) 


(3-33) 


where  9  is  the  geomagnetic  colatitude  of  the  point,  fi0  is  the  geographic 
colatitude  of  the  point,  To  is  the  geographic  longitude  of  the  point,  9^ 
is  the  geographic  colatitude  of  the  dipole  axis  ,  and  t0  is  the  geographic 
longitude  of  the  dipole  axis.  First,  rotate  the  coordinate  system  clock¬ 
wise  about  the  f  axis  through  (tt/2-DIP)  to  align  the  a-axis  with  the  Z 
(vertical)  axis-  then,  rotate  the  coordinate  system  counterclockwise  about 
the  Z-axis  through  (tt/2+DEC)  to  align  the  X  and  Y-axes  with  east  and 
north,  respectively.  Therefore,  to  transform  the  components  of  the  ion 
velocity  from  field  aligned  coordinates  to  local  tangent  plane  coordi¬ 
nates: 


V  =  -V.s inDIPs inDEC  +  V  cosDEC  -  V  cosDIPsinDEC  [cm  sec-1]  (3-34) 
x  a  $  a 


V  =  -V0sinniPcosDEC  -  V^sinDEC  -  V^cosDIPcosDEC  [cm  sec'1]  (3-35) 
Vz  =  -V^cosDIP  +  V^sinOIP  [cm  sec-1]  (3-36) 


and  similarly  for  the  components  of  the  magnetic  field  vector. 


Next,  we  must  transform  the  components  of  the  mean  electron 
density  gradient  and  the  second  derivative  from  the  curvilinear  dipolar 


coordinates  to  the  tanqent  plane  coordinate  system.  The  gradient 
transforms  as  a  vector,  and  the  tr ansformat ion  is  straightforward.  First, 
we  compute  the  linear  components  of  the  gradient  along  the  Cartesian  axes 
aligned  along  the  t,  R,  and  d>  directions: 


of 

1 

Of 

TU 

h 

la 

X 

rx 

if 

=  1 

If 

h 

3 

>f 

=  1 

If 

h 

3-i 

h  , 
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h  are 

■t> 
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x,  f 

,  and  $ 

r7_: 

h 

/ 

l 

'  1  +  3  < 

(3-37) 


(3-38) 


(3-39) 


[cm] 


(3-40) 


m 


h  ^  =  ?R. 


h  -  7s in^ 


o  +  cos  ^  -  1 


1  +  3cos?n 


[cm] 

r  cm] 


(3-41) 


(3-42) 


where  ?-  is  the  radial  distance  to  the  point  [cm]  from  the  center  of  the 
earth,  is  the  ratio  of  that  distance  to  the  radius  of  the  earth,  and  n 

is  the  magnetic  colatitude  of  the  point  from  Equation  3-14.  The  compo¬ 
nents  of  the  gradient  are  then  transformed  to  the  tanoent  plane  system  by 


102 


the  two  rotations  defined  in  Equations  3-34,  3-36.  Finally,  the  gradient 
of  the  log  density  is  converted  to  the  mean  density  gradient  using 


qrad(Ne)  =  Negrad(f)  =  N^grad f 1 n ( N  )  1  [cm-4]  (3-43) 


where  N  is  the  interpolated  local  mean  electron  density. 


The  second  derivative  is  more  complicated  because  it  is  a  second 
rank  tensor  and  must  be  transformed  accordingly.  We  must  first  select  a 
fixed  coordinate  system  in  which  we  wish  to  calculate  the  components.  The 
local  tangent  plane  is  not  satisfactory  because  its  orientation  changes 
with  a  change  in  location.  The  LOS  aligned  system  would  seem  to  be  the 
most  loqical  choice  because  it  is  fixed  and  the  components  must  ultimately 
be  transformed  to  that  system  anyway.  It  is  preferable,  however,  to 
define  a  Cartesian  system  with  the  z-axis  aligned  along  the  dipole  axis, 
and  the  y-axis  defined  by  the  cross  product  of  the  rotation  axis  and  the 


dipole  axis. 

The  coordinates  of  a  point 

in  that  system  are  just 

X 

=  Rs  inn  cos'!) 
m 

[cm] 

(3-44) 

y 

=  Rsinh  sinfi 
m 

[cm] 

(3-45) 

z 

=  Rcosn 

m 

[cm] 

(3-46) 

where  R  is  the  radial  distance  [cm]  to  the  point  from  the  center  of  the 
earth,  is  the  magnetic  colatitude,  and  is  the  magnetic  longitude. 
For  a  general  transformat  ion  from  curvilinear  to  Cartesian  coordinates, 
the  components  of  the  gradient  can  be  shown  to  be 


}f  _  r  1  ^Xn  3f 

tt;  j  h]  ^ 


(3-47) 
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3y  applying  this  operator  to  itself,  we  obtain  a  general  expression  for 
the  components  of  the  second  derivative: 


32  f 


3x  3x^ 
m  n 


3 

l 

j 

3 
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3  ,  3x  3x  , 

y  1  m  n  3f  3  ^  1  ^ 

k  h2.  3uj  3 u j  h2 


32f 


3  ,  .  3x  5x 

y  1  1  m  n 

k  h^  h2  3uj  3u)< 
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,  -j  j  ,  ,  3x  32x 

+  1  y  y  1  1  3  f  j-  m  n 

2  j  k  h2  h2  3u,  8uj  3 ui 3 u k 


3x  32  x 
n  m 


3u  .  3 u  .  3  u, 
J  J  k 


]  (3-48) 


The  derivation  of  the  coordinate  differentials  3xn/3Uj  is  complicated 
by  the  fact  that  R  and  (and  hence  x,  y,  and  z)  cannot  be  expressed 
analytically  as  functions  of  a  and  8.  That  is,  the  inverses  of  Equations 
3-16  and  3-17  do  not  exist  in  analytic  form.  We  can,  however,  derive 
analytic  expressions  for  the  differentials  if  we  introduce  a  third  vari¬ 
able.  For  a  given  a  and  8,  we  may  obtain  the  corresponding  values  for  R 
and  from  the  root  of  the  quartic  equation 


<2d;l+  +  C-l  =  0  with  0  <C<1 

where 


then 


and 


cos  a 

K  =  - 

sin48 


R 


R 

e 

sin26 


cose  =  <c2 

m 


(3-49) 


(3-50) 


(3-51) 


(3-52) 
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Equation  3-49  can  be  solved  easily  by  a  simple  iterative  techni¬ 
que.  PRPSIM  uses  Newton's  method  w’th  6  iterations.  For  k  _<  10,  this  is 
accurate  to  within  0.13%. 

Although  none  of  the  parameters  in  question  can  be  represented 
as  analytic  functions  of  a  and  ft,  they  can  all  be  represented  as  the 
products  of  analytic  functions  of  a,  ft,  and  5.  Let  us  represent  the 
unknown  function  g  as 

q(a , ft  )  =  u(a)  v(fi)  w (rj  (3-53) 


Then,  because  '  is  a  function  of  <  alone,  we  may  write 

=  V(fi)w(rj  p-  *  u( a ) v( ft )  (3-54) 

3a  da  dj;  dr  3a 

li  -  u(a)w (,)  p  +  U(a)v(ft)  (3-55) 

3ft  dft  d;;  dr  3ft 

where 

H  =  -  S1_n  1  =  -  <  tana  (3-56) 

3a  sin4ft 


—  =  -4r  cotft 

3ft 

and  from  Equation  3-49  it  is  easily  shown  that 

dc  =  -  /i  I  ^ 
dr  4-3g 


(3-57) 


(3-58) 


By  employing  this  strategy,  we  obtain  the  expressions  for  the 
derivatives  of  the  scale  factors: 


3  {  h-?,  =  _  2  +  6c2)  +  1 

3a  ^  (4-3g)2 


COSa  Sina 


h-‘ 


(3-59) 
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-  1}  h'2 
1  a 


( 3-60) 


r 


i-  (h-2)  =  -  4cot8  1 2t8~.15l;  *  60 


aR  a 


(4-30' 


3 

34 

(h-2) 
v  a  1 

=  0 

(3-61) 

3 

f  h"2 ) 

12  tana  !. 

h-2 

(3-62) 

3a 

B 

(4-3n)2 

3 

fh-2l 

2  { ?4  li 

-  r 

,(2-‘>co.b  -  ,lWs)  I  h- 

2  (3-63) 

3  8 

B 

(4 

.3^)2  cosB  sing  6 

3 

(h-2) 

-  0 

(3-64) 

3$ 
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(h-2) 

_ 

6tana  | 

(l 

JlL}  h'2 

(3-65) 
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(2 

-3c)|  h-2 

(3-66) 

33 

4> 

(4 

-3-; )  * 

3 

rK2) 

=  0 

(3-67) 

For  the  coordinate  differentials,  we  obtain 
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(4-3c) 
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3x 
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Tr 

(4-3c) 
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3x 

a* 


-  y 


H  -  X 

3  4> 


21  =  o 

3t> 


And  for  the  second  order  coordinate  di f ferent i al s : 


1  32x 
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(3-72) 

(3-73) 

(3-74) 

(3-75) 

(3-76) 

(3-77) 

(3-78) 

(3-79) 

(3-80) 

(3-81) 

(3-82) 
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32  x 
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(3-83) 
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(3-84) 
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3  a3  4> 

32x 

3y 

(3-86) 

3  S  3  4> 

38 

32y  = 

3X 

(3-8/) 

3  S  3  4> 

Is 

92z 

0 

(3-88) 

3  63  4> 

The  second  order 

coordinate  di ff erent i als  are,  of  course, 

symmetric  with 

respect  to  the  order  of  differentiation,  i.e.. 

2 

d  X 

32x 

n 

n 

(3-89) 

3u  .  3u. 

3u.  3u  . 

J  k  k  j 

The  quantities  defined  in  Equations  3-59  -  3-88  are  then  used  in 
Equation  3-48  to  calculate  the  components  of  the  second  derivative  in  the 
dipole  aliqned  Cartesian  system. 

Next,  we  must  rotate  the  second  derivative  into  the  tangent 

plane  coordinate  system.  This  is  accomplished  by  a  series  of  four  rota¬ 
tions.  The  first  two  rotations  rotate  the  coordinate  axes  from  a  dipole 
aliqned  system  to  a  rotation  axis  aliqned  system.  If  the  dipole  axis 

points  toward  colatitude  9q  and  east  longitude  rotate  the 

system  clockwise  about  the  y-axis  throuqh  0^  to  align  the  z-axis  with 

the  rotation  axis.  Next,  rotate  the  system  clockwise  about  the  z-axis 
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throuqh  ^  to  aliqn  the  x-axis  with  the  Greenwich  meridian.  The  next 
two  rotations  rotate  the  qeocentric  Cartesian  axes  into  the  tanqent  plane 
system.  If  the  qeoqraphic  colatitude  and  longitude  of  the  site  are  0O  and 
>3,  respectively,  rotate  the  system  counterclockwise  about  the  z-axis 
through  so  that  the  x-axis  points  to  the  riqht,  and  the  y-axis  is 
directed  away  from  the  observer.  Next  rotate  the  system  counterclockwise 
about  the  x-axis  throuqh  -30  to  aliqn  the  z-axis  with  the  vertical. 


For  a  clockwise  rotation  of  0  about  the  x-axis: 
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For  a  clockwise  rotation  of  0  about  the  y-axis: 

J  ^  =  cos  2  3  — -  +  ?cos  Osin  0  J-JL  +  sin  2e  JLf 
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And  for  a  clockwise  rotation  of  e  about  the  z-axis: 
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— - -  =  cos0sin0  (Ll  -  Ll)  +  (cos20  -  sin20)  I— (3-105) 
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Finally,  the  components  of  the  second  derivative  of  the  loq 
density  are  converted  to  the  components  of  the  second  derivative  of  the 
mean  density  itself: 
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